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AN AUTOMATED SYSTEM FOR INTERPRETING GRAVITATIONAL 
ANOMALIES (THE MINIMIZATION METHOD) 


Ye. G. Bulakh 


Introduction 


Geophysical survey methods are of primary importance in geological surveying. /3* 
A properly designed survey using geophysical methods makes it possible to direct 
future survey work efficiently and to decrease its cost significantly. 


One of the most important components of any geophysical survey method is 
the geological interpretation of field material. In the end result this inter- 
pretation, together with correct methodology and field operational technology, 
determines the success of survey operations carried out using geophysical methods. 


Computers have recently begun to be used to process and interpret geophysical 
materials. They have been used widely to interpret gravitational data. 


While not setting ourselves the goal of reviewing the possibilities of 
applying computers to all areas of exploratory geophysics, we will examine 
several questions concerning the use of computers to interpret gravitational 
anomalies. 


It is possible to distinguish four independent trends among questions 
concerning the possibility of using and applying computers to quantatively 
interpret gravitational anomalies. 


To the first trend we assign work associated with transforming observed 
anomalies and with calculating various correction factors. This work is usually 
very unwieldy. The results of the conversions and transformations of the field 
observed can be used both for qualitative and for quantatitive interpretation. 


The second trend includes the problem of directly interpreting anomalies 
with existing computers. This involves programming and machine solution of a 
number of tasks relating to calculation of parameters which determine perturbing 
masses. In this process the initial data are represented by the values of the 
anomalous function. 


The third trend includes questions concerning use of computers to calculate 
various measuring grids, nomograms, tables, and other interpretative aids. These 
aids can increase the effectiveness of methods of interpretation already known, [4 
and can aid in developing new methods which could not be introduced earlier 


because of the cumbersome calculations involved. 


*Numbers in the margin indicate pagination in the foreign text. 


The fourth trend comprises work toward creation of specialized problem- 
solving machines to interpret anomalies. 


Let us examine each trend in somewhat more detail. 


First trend. Computers make it possible to perform calculations relatively 
quickly even with the most complicated formulas. Hence it is natural that, from 
the very start, research on the use of computers to process and interpret geo- 
physical observations was aimed at automating laborious calculations. 


By the time computer technology was introduced into gravimetry, the 
calculations associated with field transformation were extremely cumbersome. 
Therefore, programs were written primarily on field transformation to perform 
gravimetric calculation problems. 


In the USSR I. A. Balabushevich, B. V. Bondarenko, R. S. Volodarskiy, 
O. K. Litvinenko, V. N. Strakhov, and others became the first to use computers 
for the transformation of fields. The use of computers for the transformation 
of anomalies required the development of more sophisticated computer design 
(G. I. Karatayev, A. K. Malovichko, M. G. Serbulenko, V. N. Strakhov, M.-LaPorte, 
and others). This trend includes work in computer calculation of correction 
factors for the influence of landscape relief in surveying with a gravimeter 
(V. I. Aronov, L. A. Koval', M. Bott, M. Kane, and others). A fairly circumstan- 
tial review of the research in this direction has been made in [36, 49, 85, 111, 
etc.]. 


The second trend, i.e., use of computers to interpret anomalies, was developed 
almost simultaneously with the first trend. Computers make it possible to carry 
out even the most difficult calculations in a relatively short time. However, 
the basic importance of these investigations for the interpretation of anomalies 
is the fact that, with the new capabilities of technology, qualitatively 
improved methods of interpretation can be created. Now there arises the task 
of developing more universal methods of interpretation which cannot be applied 
without using computers. It is impossible to create one universal machine 
method for interpreting gravitational anomalies; the interpretation methods 
differ under various physicogeological conditions. It is necessary to generate 
a small array of these methods and then carefully determine the areas in which 
they can be used. Works are already available on creating effective machine 
methods of interpretation. 


The first stage of the second trend is development of effective methods 
of solving direct problems, a task which has great significance for automated 
computation in interpreting anomalies. No matter the methods by which anomalies 
are interpreted, the diagram of the geological structure of a region will be 
built, as a result, on the basis of all the data (geological and geophysical). 
It is necessary to solve the direct problem of gravitational survey to confirm 
the correctness of this diagram's structure. In addition, many researchers are 
coming to the conclusion that when using field transformation it is impossible 
to single out an anomaly caused only by the geological feature in question. 
The so-called geological reduction method has been widely applied in the 


practice of interpretation. Its basic characteristic is the fact that the 
fields are divided by subtracting from the observed gravitational field the 
effict of the feature, the position and dimensions of which are established 
according to geological data or during interpretation of the data of other 
geophysical methods. In this process it is necessary again to solve the direct 
geological structure problems formulated by the researcher. 


When calculating anomalies for certain geological features, we encounter 
an abundance of perturbing masses. This causes great difficulties. However, 
the calculation process for solving this direct problem can be unified with an 
adequate degree of precision. 


In fact, say that it has been established to be possible to consider a 
perturbing body as being two-dimensional. Let us further assume that the 
perturbing masses are concentrated in one or several areas. Each area can be 
well approximated by a segmented straight line boundary. Thus, the gravitational 
field created by any two-dimensional body is approximated by the sum of the 
gravitational fields of its projections. Each projection is determined by five 
parameters: X,> h, Xo» H, which are the coordinates of the angular points, and 


o which is the excess density. 


Calculations for three-dimensional bodies can be unified in a similar way. 
For this purpose it is sufficient to represent the perturbing bodies as the 
sum of projections of limited extent or as the sum of parallelepipeds. In the 
latter case each geological feature is described by seven parameters. Six of 
them characterize the location of the boundaries, and the seventh characterizes 
the excess density. This approximation makes it possible to describe bodies 
of variable density. 


The programs which have been written make it possible to calculate the 
observed field. From the interpretation stage, in which the hypothesis 
concerning the geological structure of the region was formulated, one proceeds 
to calculate the anomaly. It is necessary to feed into the computer information 
which reflects the determined geological configuration. This is a system of 
five-dimensional vectors (for a two-dimensional geological section) or seven- 
dimensional ones (for a three-dimensional section). In addition, it is necessary 
to have a clear idea of the coordinates of the initial point, the distance 
between the points on the profile, and the distance between the profiles, as 
well as the coordinates of the terminal point. The programs are written for 
an indeterminate number of vectors which reflect the geological structure. A 
maximum limit is imposed on this number, depending on the volume of the effective 
memory. Thus, it is necessary in each specific case to feed into the computer /6 
a number equal to the number of vectors providing information which describes 
the geological configuration. 


The effectiveness of these calculations has been verified by means of 
numerous examples in interpreting regions complex in the geological sense of 
the word. 


O. K. Litvinenko, M. Bott, M. Tal'vani, and others have concerned them- 
selves with the solution of similar problems. 


The second stage is the development of automatic machine methods for 
interpreting anomalies. These methods have been given a theoretical basis by 
the research of A. N. Tikhonov, V. K. Ivanov, M. M. Lavrent'ev, and others. It 
is also possible to classify the methods of searching for singular points [59, 
60, 74] under the works in this direction. Relatively stable algorithms have 
been obtained which make it possible to carry out the calculation processes by 
computer. 


In the works of S. V. Shalayev [75-77], study is made of the problem of 
approximating an observed anomaly by a rational fraction. Such an approximation 
makes it possible to find the singular points of perturbing masses. A linear 
programming device is used for the actual calculations. 


In [56a, 80] a summary is given of the methods of determining the position 
and dimensions of perturbing bodies and the excess mass (the grid method). A 
Similar task situation is also presented in [89]. The half-space in which the 
perturbing masses are concentrated is divided into elementary bodies (prisms 
or parallelepipeds). By minimizing the difference between the observed and 
theoretical anomalies, one finds the values of the excess densities within each 
elementary body. S. V. Shalayev [77, 78] developed this method further in his 
works. He reduced the problem to linear programming. 


A second variation of the grid method has been summarized in [31]. A 
number of works have since appeared, in which this variation has been put to 
practical use [37, 52, 53]. 


In a number of works the solution of the inverse problem is reduced to 
iteration processes [87, 91, 114, etc.]. In order to minimize the specially 
constructed functional, the least square method [90, 93, 96] is used. 


In [100] the relaxation method is used to solve the inverse problem. The 
perturbing body is approximated by the sum of vertical prisms with their square 
bases. The gravitational effect is calculated at points of a square grid. The 
method can be used if one knows the depth of the interface at least one point 
and the excess density. 


On the basis of the Poincaire theorem, D. Zidarov has worked out an original 
method for determining the configurations of perturbing geological features 
[117]. If the centers of the masses are known, then it is possible to construct 
a convergent iteration process for "sweeping out'' the perturbing masses into 
the surrounding space. The computation cycles are repeated until homogenous /7 
geological features are obtained which have the assigned excess density. 


Methods of special correlation analysis [33] have been developed for fore- 
casting in geology. By taking a certain area as a standard, it is possible to 
predict geological structure on the basis of observations of geophysical fields 
in other similar regions. 


I, Nedyalkov [102] has proposed a heuristic method for solving one class 
of inverse problems of the theory of potential. The parameters of the geological 
configuration are determined by using specially constructed self-learning [sic] 
algorithms. 


One could continue this list of the works the theme of which adheres to 
this trend. They would include the fundamental monograph by F. M. Gol'tsman 
[20] and a large number of works by other researchers. The works cited in this 
review can serve only to illustrate the general content of the second trend 
in research. 


As regards the third trend, as indicated earlier, it comprises works which 
make use of computers to calculate tables, nomograms, measuring grids, etc. 
The results of the calculations serve the purpose of practical application of 
a new method or a new means of processing and interpreting anomalies. At first 
glance the question of the means used to calculate these measuring grids or 
tables is not one of fundamental importance. This is true when the formulas 
are not too cumbersome and allow calculations using tables and an adding 
machine. However, one often encounters problems solution of which cannot be 
accomplished without using computer technology. The calculation of integral 
measuring grids [14] may serve as an example. 


It is not possible to give a complete survey of works in the third trend 
of research. Many calculations have been carried out recently using computers, 
and by no means all authors concentrate their attention in this area [82-84]. 


A review of the works in the fourth trend has been given in [57]. 
Almost all of the specialized analog devices are based on the fact that the 
various physical fields have common properties. These devices make use of the 
Mathematical analogy between gravitational fields and electrical, radioactive, 
electromagnetic, and light fields. 


Researchers have come to the conclusion that analog computers can be used 
effectively to interpret an observed field only when the following conditions 
are satisfied: 

1. The parameters of the perturbing body undergo slight variation; 


2. Calculation of the field is done relatively quickly; 


3. The results of the calculation are printed out in a form in which the 
calculated and observed anomalies can be quickly and easily compared; 


4. The device must be simple, reliable, and convenient to use. 


Despite their seeming simplicity, the models set up to interpret by the /8 
matching method have not been widely used. One of the basic reasons for this 
is the fact that by no means all of the analog computers meet the above- 
mentioned requirements. 


The problems in the second area of research have been examined in this 
monograph. Research is cited which uses the minimization method to solve the 
inverse problem. This method makes it possible to create a closed automated 
system for quantitative calculations. This system of calculations includes the 
stages of the preliminary survey type (calculation of the direct problem), 
search for a background effect, and lastly automated search for a new variation 
of the geological configuration. 


CHAPTER I 


SOLVING INVERSE PRELIMINARY GRAVIMETER SURVEY 
PROBLEMS BY THE MINIMALIZATION METHOD 


Section 1. Statement of the Problem /9 


When we speak about using computers to process and interpret gravimetric 
data, we mean primarily the automation of calculation operations. Computers 
unquestionably make it possible to carry out even the most complicated cal- 
culations relatively quickly and accurately. The basic reason for using these 
machines in interpretation practice is nevertheless the fact that, on the basis 
of new technical capabilities, quantitatively new methods of processing and 
interpretation can be developed and used. 


The entire cycle of calculation operations and logical operations can be 
reduced to a single automated system for processing and interpreting gravimetric 
data. Structurally this system consists of two independent parts. The first 
part consists of the processing of observations, the introduction of various 
correction factors and reductions, and the process of obtaining an anomalous 
gravitational field. If necessary, various transformations can be carried out 
here. The field can be converted to any other level, higher derivatives or 
field potentials can be computed, and other special functions can be constructed 
for the field under study, such as the Saksov function, the Nigard function, 
etc. The results of the processing are usually presented in the form of 
catalogs and graphs along fixed profiles and charts of the anamalous field and 
its transformations. After comprehensive processing of the observed data, 
it is possible to thoroughly analyze all the material obtained and compare it 
with the data of other geophysical methods and with information on the geological 
structure of the region. This is the total of the subject of qualitative in- 
terpretation. 


Various reductions are introduced into the observed field; correction 
factors or observed data have been converted to other functions so that it will 
be possible to state with confidence that the anomaly obtained has been caused 
only by heterogeneities in the geological structure. On a plan or graphs places 
are pinpointed in which bodies with too little or too much mass are concentrated, 
the configurations of these masses are established in general form, their 
approximate epicenters are indicated, and several conclusions are drawn about 
the geometric form of the geological features. 


The second step is represented by quantitative calculations. The numerical /10 
values of the parameters of the geological features and the density characteristics 
of the rocks which make up these features must be established on the basis of 
the peculiarities of the anomalies. The choice of the methods of making the 
quantitative calculations depends on both the type of anomalous field and on 
the geological premises which the researcher has accepted as his basic working 


material. 
7 


The second part of the automated system, i.e., the system for interpreting 
the gravimetric data, is associated with the quantitative interpretation. The 
magnitudes which characterize the occurrence of the masses of interest to us can 
be estimated quantitatively by several methods. Each of these methods has 
both strong and weak aspects, and as a result of this each one can give satis- 
factory and unsatisfactory results in determining the indicated magnitudes. 

This depends on the geological conditions under which the anomaly being inter- 
preted has been obtained. Where one of the means gives satisfactory results and 
can therefore be used successfully, another one is not used. 


As a whole all of these methods should be considered to be of equal value, 
but each of them must be used while taking into account the geological structure 
of the region in which the anomaly being interpreted is located. If the region 
being studied is relatively complicated in the physico-geological sense, and 
if the researcher has somewhat limited information about the structure and 
makeup of the geological objects, then often the only effective means of in- 
terpreting under these conditions is the selection method. The researcher 
constructs his geological model on the basis of all information available about 
the region's geological structure, while taking into account the anomalous 
field. The direct problem is solved, and the gravitational effect is found 
from the appropriate geological model. By comparing the observed and theoretically 
calculated anomalies, the researcher should change the geological model con- 
structed earlier in such a way that the newly calculated gravitational effect 
comes as close as possible to the actually observed field. In this process it 
is first necessary to know how to solve the direct problem quickly, and second 
to know how to construct the type of algorithm which would make it possible 
to answer the questions of how to use the geological model so that the observed 
anomaly coincides in the best way possible with one which has been calculated 
theoretically. 


In the following section we will dwell on a description of this type of 
algorithm. Calculation of the direct problem is a component part of this 
algorithm. Thus, in order to solve the problem of determining the parameters 
of geological bodies, we have an observed anomaly, all variations in which are 
related to geological heterogeneities and the initial variation of the geological 
structure model. The observed field and the tentative model of the geological 
structure together form the initial data for the quantitative calculation. It 
is necessary to change some of the geological parameters so the observed and 
theoretically calculated fields may more closely coincide. 


Now in the observed field let us fix end points with the coordinates (x; /1l 


and y;)- Let these be the most widely varied points at which the observed 


gravitational anomaly manifests itself in its most characteristic manner 
(extremums, points of inflection, etc.). In the future we will select a field 
only at the fixed point. 


We will introduce one restriction. Let the geological model be constructed 
in such a way that its gravitational effect can be recorded in an analytic form. 


If the geological bodies are approximated by putting together projections (gra- 
vitational steps) which are finite in space, then even the most complicated 
structural model can be assembled with relative ease from these elementary 
features. 


Thus on the one hand we have the observed anomaly V (x, y,), and on 


obse 


the other hand the theoretical anomaly Viheo Py? Poo+++Prs X55 y;,)> 1S); 2, 


..,n. Here the function V can stand for any component of the gravitational 
field (Ag, V xz? etc: ). If the geological model has been described using m 


elementary bodies, and each body is characterized by s parameters, then the 
result of solving the direct pechien welt be the theoretical anomaly 


ve Gs y= SVP, x,y). (1.1) 


CO ao f 


ame a 


By the symbol Po we indicate the vector which characterizes the location and 


dimensions of the elementary bodies, P, = (P5y> Pigs -++Pyg)> Let us illustrate 


this with an example. Let the geological structure of an area be approximated 
by a collection of direct projections which are limited in space. The location 
and dimensions of each projection are characterized by these parameters: o is 


the excess density, h and H are the depths to the upper and lower limits, Z7 and 


Lo are the parameters which characterize the projection's dimensions in space, 


and d is the location of the vertical boundary relative to the selected beginning 


of the coordinates. The value of these parameters is shown in Figure 1. Thus 
the totality of the vector P, = (G53 hig He, 2.25. byes d,) characterize the 


Jo ed 29 
geological schematic of the region being investigated. 
Now it is necessary to find the value of the parameters p, at which the 


J 
function V best coincides with V The process of comparing these two 


theo obse’ 
curves can vary. Let us examine one of them. 


In order to compare the functions Vobse & y) and Veneol® y), let us set 


up the functional - 
ERE Vs. Vobealie 4) — Vico 9). (1.2) 


Since the points with the coordinate (x, Y;) have been fixed, then F depends 
only on the parameters p,. Having numbered them sequentially, one can write 
J 
TF =F (Pty Paw +--+ + PN) ga) 


Expression (1.3) can be regarded as a function from one N-dimensional vector 
P= (Py> Poo eee> Py) - 


/12 


Let us select a vector P so that the 
functional (1.2) acquires minimum value. 
Thus the task of interpretation is reduced 
to the extreme problem of the function with 
many variables. 


Let us turn our attention to some 
theorems about the extremums of several 
variables. The following theorem speaks of 
the condition necessary for the existence 
of an extremum. Let function (1.3) be such 
that in the area where it is being investi- 
gated there are partial derivatives of the 
first and second orders. Then if at point 

(o) _{0) 
M, (Py °> P5 
exists, then all partial derivatives of the 
first order at this point are equal to 0. 


, lll, pi?) a local extremum 


Figure 1. Values of Para- 


meters Characterizing the Then 

Location and Dimensions of el ee 

an Elementary Geological Body, “Tear li” (i=l, 2,..., ) (1.4) 
i.e., the Direct Projection HS 

Limited in Space. is the necessary condition for the local 


extremum. System (1.4) makes it possible 
to find the points of the possible extremum 
(or a Stationary point). It is possible to have cases where the function 
F (Py> Pos sess Py) does not attain the extremum, and thus (1.4) turns out to 


be a necessary condition, but not a sufficient one. 
In order to obtain an unambiguous answer to the question of extremums, it 


is necessary to use another theorem. In it another differential of function 
(1.3) is examined. This differential is the quadratic form from the differential 


dp, dp., ee dp, and can be written in symbolic form as 
Lee oe a ay eee 8 2 
PF = (ape det ap aPat «+ + Gp dw) F (1.5) 


and in general form as 
oa ~ N 
OF = O= 2, ayhAp ay = Aji- (1.5a) 


Here the differential dp, is designated by h The quadratic form (1.5a) is 
Py gn q 


k' 
called positively determined (negatively determined) if, for any values of the 
variables not equal to 0, at the same time this form has positive (negative) 
value. The positively and negatively determined forms are combined under a 
general designation, i.e., fixed forms. If the quadratic form (1.5a) acquires 
both positive and negative values, then it is called sign variable. If in 

(1.5a) the form has values of one sign, but at some point acquires 0 values, then 
this quadratic form is called quasi-fixed. 


10 
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Let us now turn to the theorem itself. Let the extremum of function 
F =F (Py> Pos sees Py) be possible at point M, @,, pe atone ae Let 
us examine at this point the second diceseential 2°F/m, . This differential can 
be written ae expression (1.5). If (1.5) at point M, has a local extremm. 
= aes if a? F < 0, then the function at point M Pe a Maximum, and if 


d oe > 0, then it reaches a minimum. When the second differential is a sign- 
variable form, then at this point the function does not have an extremum. When 
the second differential has a quasi-fixed form, the function can either have or 
lack an extremum. Further research is required to clarify this question. 


The criterion for determining the sign of a quadratic form was established 
by Sil'vestr and is named after him. 


The quadratic form (1.5a) is characterized by the symmetrical matrix 


(1.6) 


and are called the principal minors of matrix A. 


In order for the quadratic form (1.5a) to be positively determined, it is 
necessary and sufficient to fulfill the following inequality 


ee > 0, Ay >0, yAy>9, . , Ay 0. - 


In order for the quadratic form to be negatively determined, it is necessary and 
sufficient that the signs of the principal minors Ay, A » ., A,, alternate 


N 
while Ay < 0. 
This is the classic method. To implement it, it is necessary to find a 
solution for system (1.4). When solving gravimetric problems, this is a /14 


system of transcendental equations. This system is not solved in a general form. 
In every actual case it is solved by using approximate numerical methods. 


It is possible to approach the minimization of functional (1.3) in another 
way. Some restrictions in the form of 
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Or Pr Po - ss ’ pu) <0; t= 1, 2, eee, A (1 ; 7) 


can be imposed on the parameters pss The area D is determined by the system of 


relationship (1.7). The problem consists of finding among the points of area D 
the point P* for which the following is attained © 

F(P*) = minF(P), P€D. 
The problem formulated in this way relates to the class of problems in non- 
linear programming. In its most general form this is a problem of functional 
programming. Depending on the type of function F and the area D which is 
determined by the inequalities (1.7), the problem can be simplified by reduction 
to problems for which solution methods are known. In particular if the function 


(1.2) is convex and smooth, and if the inequalities (1.7) determine the convex 
area D, then this problem is reduced to convex programming. 


The concept that the function is convex is very important. For convex 
functions, theorems have been established regarding the criteria for the 
uniqueness of the inverse problem. Let us examine some formulas and theorems 
of convex function. Here they have been presented without proof. A more 
complete analysis of this question can be found in [13, 27]. 


Let two points pt) and p (2) be assigned in N-dimensional space. The set 


of points P = p(1) + cp?) - pl), where t is any number within the segment 


[0, 1] (0 St <1), is called section p() p(?) connecting these points. 
(1) p (2) 


Set P is called convex if, together with any two of its points P 


> 


all points in the section p(t) | p(?) belong to this set. 


Let the function F (P), P = [p,> Poo ses Py! in N dimensional space. 


F (P) is a convex function in set P if for any two points p(1) and p (2) the 
following condition is satisfied 


FP! 4 1(Po—P)) <F (PY) + (1.8) 
op ELE (POS F (PM): 


If in (1.8) the inequality sign is strictly correct, then one speaks of a strict 
convexity. 


The theorem which established the criterion of convexity is extremely 
important. Let the function F (P), which is assigned in the convex set P be 
differentiated twice. If at two points in this set the second differential 


De. 3 eee : : 3 F 
dF is a positively determined form, then F (P) is a convex function in set P. 
Thus on the basis of an investigation of the second differential, one concludes 
that the function is convex. 


Why is it so important to establish the convexity of the minimized function? 
The reason is the uniqueness of the minimum. The following theorem confirms this. 
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The differentiated, strictly convex function F (P) which is assigned in 
the convex set P can have a local minimum only at one point in this set. 


It is true that, in its most general form, functional (1.3) can be non- 
convex since the inverse problems fall within the class of incorrect problems 
and do not have a unique solution. However, by imposing restrictions on the 
form of the elementary body and by assigning some of its parameters, we can 
always ensure that the solution to the problem will not come from a particular 
class. This unique controlling factor ensures that the sole solution to the 
problem has been achieved. 


Two basic tendencies are distinguishable among the methods for solving 
non-linear programming problems. One of them combines the methods of a systematic 
search for a solution. Various multi-step searches for a solution are constructed. 
The direction of the search is selected at each step, while the various local 
properties of the minimized function (for instance, the gradient of this func- 
tion) are used. The second tendency uses the random search method. Each 
tendency has both strong and weak points. 


The systematic method has complex algorithms. If the minimized function 
has local or boundary minimums, then similarity is, as a rule, not guaranteed 
between the results of the calculation and the global principal minimum. A 
sampling of the initial points of the search is required in a search for the 
global minimun. 


In the random search method the algorithms are relatively simple, but they 
require that the entire function be calculated frequently. 


When minimizing (1.2) or (1.3), their relative cumbersomeness should be 
kept in mind. Derivatives of the functions are computed almost without further 
expenditures of time, i.e., when computing the function itself. However, not 
all of these distinctive characteristics are germane to choosing a minimization 
method. One should take into account the fact that the choice of an initial 
approximation made by a qualified specialist who takes into account that which 
is already known about the geological peculiarities of the region under investi- 
gation makes it possible to obtain a solution to the problem in a previously 
selected class. This latter is the decisive factor in choosing a minimization 
method. In order to minimize (1.2) or (1.3), we use the gradient method of the 
fast descent. Let the vector P = (Py> Pos sees Py) turn (1.2) or (1.3) into 


a minimum. Then at all other proximate points which are characterized by the 


(k) 


vectors p~ °, the functional 


P(PM)> F (P). 
the search for the point P comprises the following: Let the initial approximation 
p (0) = Ger, ee seek Bee be known. Let us construct the vector grad /16 


° 
F cp! dy = (Fo its Fon) It is well known that the function F (P) will in- 


crease in the direction of the vector gradient, i.e., in the opposite direction 
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this function will decrease. Through point p() in the direction of the vector 


gradient, we will draw the straight line 

ae ea 

(Pa BO herad F(P), (1.9) 
By ascribing the various values to the coefficient A, we will obtain various 
points on the straight line (Figure 2). If 4 is greater than 0, then the point 
P will be located in the area where the function F (P) decreases. The value 
of the minimized function at ay point on the given straight line can be written 
thus : =, 


7 - rp = ames 


or 
We will find the value A = a at which F OC) = min. For this it is 


necessary to solve the equation 
“hon (1.10) 


Thus the point on straight line (1. 9) has been determined 
“POs P 2, grad F(P™), 
at which F cpt )y assumes the least of its possible values. Now we take the 


point p(l ) as the initial approximation and determine the new point pl?) All 
of the calculations are reduced to the following determination of the vector 
components. The calculations are carried out according to the formulas 


ptth ox pi — An (Fo) (1.11) 
ane eee ade. 


where k = 1, 2, ... - the nibtaeton wanbees First it is necessary to calculate 
the coefficient AL? If difficulties arise when setting up and solving equation 
(1.10), then it is possible to make use of the approximated value Ay as deter- 


mined by the Newton method. In contrast to the exact value, the approximated 


value of the coefficient is designated by dunt 
f Raw eee (1.12) 
FY + Fp, ie ieee Ua 
Great interest is aetpibuted to the value of the function F = Fein at 
which the minimization process can be concluded. We will assume that a mistake /17 


in calculating the function will be due only to a mistake in the observations, 
and that all calculations are carried out with considerable precision. In this 
case the error can be replaced by the differential of the function 


AF =2% [Vobse (s+ Ys) —VihealXo 9,)) AV obse 
eee, AS 
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Let us further assume that the difference between the observed and the cal- 
culated anomalies at the end of the computation will not exceed the observation 
error. 


fin — AF = = 2 p> (AV ose >. 
eR. 7 tut hig ee 
or ; a 
Re N20 (AVo bse). (1.13) 


It must be noted that all of this is correct if the perturbing geological masses 
can be described exactly using approximating bodies. In actuality the equality 
(1.1) does not approximate the observed function with absolute exactness, and 
the minimum of function (1.2) turns out to differ from 0. If the small value of 


the error in the field AV is chosen, then F_. > AF,._ can be obtained. In 
obse min fin 


this case it is necessary to minimize function (1.2) to its smallest possible 
value. One should expect that, having obtained the smooth decresing sequence 


fin. Motte FO) wil differ slightly from F‘K * 1), 


It is worthwhile to stipulate the termination of the calculation according to 
the criterion of the relative difference 


ke) pet) 
a 1 =A.: (1.14) 
yg sas 


The value A can be determined from the model investigations, or it is established 
during the process of sampling the actual practical problems. 


RD, | we will not attain F 


2. Other Variations of the Problem 

There can be other approaches to comparing the 
observed and theoretically calculated functions. 
Let us examine several of them. 


1. Let us set up the functional 


F = 3 |Vobse(* 4) — ¥gheolXi> 41) |- (1.15) 
_ he 


The latter expression can be written in a some- 


Figure 2. Search For what different manner. (1.15) includes N components, 
a Point Along a Vector each of which is taken at its absolute value. We 
Gradient. will divide these components into two groups. In 


the first group we will include those where 
- = 
LV sees (x; , Y;) ee (x, y;)] #0, and the rest 
will be in the second group. We will designate by r the number of components /18 


which have turned up in the first group, and thus the number of components in 
the second group will be (n - 2): Now (1.15) can be written 


F= > WVobsel 4.) —Kenedl te UP) | + (1.15a) 
se 


aes : We “Deneck zi Ww Vopse (4 y;))- 
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Now a specific function (1.15a) will underlie minimization. Here the 
gradient method of the fastest incline can be used without any kind of re- 
strictions. 


2. When comparing the two functions Y hee (x, y) and Moiese (x, y), we 
select one of the fixed points (x; 5 y;,)- For it we will write equation 
= Vobseltes Yi) — neo (Xs Yow P)s (1.16) 


where, just as earlier, P is the vector which characterizes the location and 
dimensions of the geological bodies P = [p,> Par sees Pyl- 


p(o) - (0) (0) 


In the initial model of the geological structure, Py os Pp os eres 
ieee One must find the values of P such that Ai is as small as possible. 
tet P = PO) + ap, ice. 


(1.17) 


Now (1.16) will depend only on (Ap), APs, Sen Apy)- As a rule the function 


Vee, (s, y, P) is relatively complex, and the parameters P; belong under the 
Signs of transcendental function. We linearize the problem. Let us put Vohes 
(P) (at a fixed point) (x, y,J) into a Taylor series and let us restrict 
ourselves to only the linear part. 
ha av vi, e 
Viheol?) # Vineo( P®) - Op, pals Ap,+ -:- — (1.18) 
} -* “ ‘4 
np. Apy. 
i a a 7 
Expression (1.16) is rewritten thus: 
-_— N 
Hiagonetes 00 4) —Vineo(P™s Xi» y)— 2 A,jAp;. (1.19) 
Here sree - 
eet av 
Ay = Ay (is y,, P° y= “Opi ‘O° 
having accepted A + 0 (1.19) can be written in the form 
BS et (1.20) 


>» A,Ap; = Voreel i ¥)— Vineo Yi» p a F 
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Parameter i acquires values equal to 1, 2, ..., nm, i.e., (1.20) is a 
system of linear equations with N unknowns. If n > N, then a redefined system 
has been obtained which can be written in matrix form as 


AAP = 6 
ne a tars ce Es Oa Sete ets? 
AnAis eae Ain Ap, \ 
A= Ay Ag wee Aon AP « Ap, . ; 
Ant An wee Aan, Apy_ 


This system can be solved by using the least-squares method. The method 
itself will relieve us of the necessity of investigating the compatibility of 
the selected system [45]. 


Without citing any calculations, we will simply indicate that vector AP is 
found from solving the system 


| A’AAP = A’, (1.21) 


where A' is the transformed matrix A. As is well known, this relationship 
always leads to a fully defined system of exactly the same number of equations 
as we have unknowns [45]. 


3. It is possible to introduce additional limitations into the function 
which should be minimized. These limitations will play the role of regulating 
parameters (according to A. N. Tikhonov). We will require that during minimiza- 
tion the desired parameters do not differ so sharply from their initial values. 
For this DUTpOSe we introduce into the functional (1.2) another member 


Sher ai where A. is the constant coefficient which plays the role 


of a weighting function. Wherever the version of the parameters, according to 
the geological data, should be small, the coefficient ae will be selected by 


the researcher to be large. If the variation of any parameter is not restricted 
by the geological data, ue can be selected to be sufficiently small or equal to 0. 


The basis of the final minimization is the functional 


a om - 


#- > = Palbobse Mn ye btheol*i gp + > A; (i — Py. (1 ° 22) 


This functional eee not only on the seeae heey parameters P; , but also on 
the regularization parameters has 


The minimization method can be used to interpret various anomalies by 
postulating differing forms for perturbing masses. Relatively simple relation- 
ships for computers calculation can be obtained in cases where the perturbing 
bodies are approximated by a group of spheres, right-angled prisms, etc. Through 
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identifying the area of the perturbing body by the total of its projections, it 
is possible to obtain a solution to the problem for an arbitrary contour with 
an arbitrary distribution of masses. 


3. Solving the Inverse Problem for a Group of Cylinders by the Anomaly es 


We will illustrate with an example the interpretation of gravitational 
anomalies using the minimization method. Let an anomaly Nae be assigned and 


let there be established the possibility of approximating the perturbing bodies 
with cylinders. The position and geometric dimensions of each cylinder can be 
characterized by the following parameters: dis the abscissa of the epicenter, 
h is the depth at which the axis occurs, and o is the excess density (Figure 3). 
Here it should be noted that the derivatives have very small values according 

to their mass. Therefore, changes in this parameter are small in comparison 

to the others. This kind of inequality in changing the parameters leads to a 
very slow convergence of the calculations' iteration processes. We will express 
the mass by means of a linear parameter. 


The mass of a unit of the eyhindey length is 
determined by the relationship M = nR*o. As is well 
known, it is not possible to determine o and R 
“= separately. We will designate oR? = t? and we will 
4 determine parameter t in the future. In order to 
be able to take into account the sign of the excess 
density, we accept 


Figure 3. Values of Se ey 
the Parameters Charac- ‘A ojRj =8 ~egre 
terizing the Location Here i eae 
of a Cylindrical Body. 


wed. . 


o- : if a>0, 
need={ i, ie get 


Thus the position and dimensions of each cylinder 


determine the following four parameters [t, h, d, sign (c)]. We will set up 
the function 


7 2 
Fan = Sfranncx— ea S| sien tay) pte: — i) I. (1.23) 


each! ad alain LE 
= [x — di + IP 


We will consider the parameters of sign (7) to be constant, and then m 
three-dimensional vectors will be unknown: 


P. = [t., h., d.]. 
J J J jl 


We find the components of these vectors, of which function (1.23) will 
become the smallest possible. 
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We will assign the initial approximation: iia le ee a”, Sign 
to) I We will determine the following approximations by’ the formulas 
| : ger) — i Aa (Ft) (1.24) 
ro & At) ni — de (Fr, ‘i 
i Lap <A (Fade Gab 2 ey me 
The derivatives will be expressed Ae f2l_ 


oh in} ~3 (x; i 
[er — i)? + AE: 


Fas = 3 a) nk vs, 


i=l 


eS Sees iy St 


where 
detergent 
8, = Vaz obsel*:) — Vaztneol%i- oe 
In the last relationships 2 Asie & 
ee thy (x a)" 
i "Ver theott) = 45 ene y sign (o;) if eae. ’ 


We will calculate the coefficient » by the Newtonian method of Newton 
(1.12). 


4. An Example of Solving the Problem 


We will examine the use of the method in the following examples. Let an 
anomaly of the horizontal gradient of the force of gravity (Figure 4) be given, 
and let it be established that the perturbing geological bodies can be related 
to two-dimensional bodies, for instance according to the method of A. A. Yun'kov 
[81]. Even a cursory analysis shows that the perturbing masses are scattered 
along the profile. It is possible to distinguish three anomalously shaped 
objects. Let us assume that the problem has been posed of evaluating individually 
each geological body, i.e., determining its mass and center of gravity. It is 
possible to approximate each geological body as a cylinder in order to solve 
the problem. The parameters of three cylinders (9 magnitudes in all) form the 
basis of the determination. Let us fix the most characteristic points of the 
observed anomaly along axis Ox. In all, 13 points are distinguished, which 
have been reduced to tabular form (Table 1). Taking into account the dimensions 
of the anomaly, meters have been adopted as linear units. Using the most general /22 
premises as a basis, we will select the model of the first approximation, the 
parameters of which have been cited in Table 2. (Some remarks of a methodological 
nature on the calculation of these parameters will be given in the second chapter). 
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| 20 eoetvoes 


The data in Tables 1 and 2 are the 
initial data for solving the problem. 
The results of the solution are given 
in Table 3. Here all the intermediate 
results of the calculations from iteration 
to iteration (some iterations have been 
omitted) are shown. 


The function at the initial values 
of the parameter was found to be equal 
to 18768 eoetvoes*. The function F in 
the following iterations acquires these 
values: 12449, 7369, 4376, 2780, 


We will determine the value of F at 
which it would be necessary to conclude 
the computation, Let us turn to formula 
(1.13). In our case n= 13. If we 
accept that the error in the observation 
is 3 eoetvoes, then F = 234 eoetvoes?. 


fin 
Figure 4. Anomaly Vw as Caused In the 26th approximation Poe = 153 
by Three Perturbing Bodies. eoetvoes“, and the value of the unknown 


vectors is Py = (147, 199, 302), P, = 


(82, 137, 600), and P. = (148, 198, 1049). After the 37th approximation, we 
obtain F = 31 eoetvoes- and Py = (149, 200, 301), P, = (86, 144, 600), and 


37 
Pz = (149, 198, 1050) 


Let us assume th 


0.8, and O25 1. In this case it is easy to calculate the radii 


We obtain Ri = 301, R 


at the cylinders have excess densities 2 te l, o 


> = 97, and Rz = 149, 


The cited anomaly was calculated for three cylinders which had the following 
parameters: Py = (150; 200; 300), P, = (89.4; 150; 600); and P. = (150; 200; 


1050). 


No. of 
the point 
rr 


20 


TABLE 1. 


200 | 300 } 400 | 500 | 550 | 600-| 650-|.800 | 950 | 1050 | J150 1350 
168 | 32 |—s8|—s| 13 | 38 fas | 159] 123 | -17 ~162) a 


TABLE 2. 
No. of the |. -, 
Be airbing t h a 
od 


J 


1 120 | 200 | 320 ° 
2 | 60 | 100/570 
3 130 | 220 |1040 


TABLE 3. 


Parameter No. of Approximation 


of the 


Exact 
value 


Perturbing | 

Body 
4 120] 123 180 
hy 200 | 198 200 
4, 320 | 319 300 
t 60 | 54 89,4 
hy 100 | 104 158 
4, 570 | 573 600 
4 130 | 134 150 
hy 220 | 218 201 | 201} 201) 201} 200} 200] 200] 199] 199] 193} 198] 200 
d, 1040 | 1040 ; 1042 1049 | 1049 | 1049 | 1049 | 1049 | 1049 | 1049 | 1049 | 1049} 1049/ 1050] 1050 
F 18768 | 12499] 7369 | 4376 | 2740 | 2378 1709 | 1092 | 909 | 1098] 737] 823] S78} 603] 396] 153] a1} — 


5. Nature of the Convergence of the Fast Descent Method When Solving Inverse 
Problems For a Group of Cylinders 


Let us examine Table 3 in which the results of the computations are cited. 
At first the function decreases. In the seventh iteration a disturbance is 
observed in the function's smooth progression. Meanwhile, the gradient of the 
function's change in the following step acquires its largest value. After the 
10th iteration, one again observes an increase in the value of F which must be 
minimized. Sudden changes in the values of the function become relatively 
frequent from this point on. 


What is the cause of the sudden changes? 


As was noted earlier, a vector gradient is sought by the fast descent method. 
The greatest change in the function occurs in the vicinity of the selected 
approximation along this vector. At the very beginning the change is sharply 
reduced, and then, reaching a minimum at some certain point, it begins to in- 
crease. The problem is to determine the point at which the function attains the 
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minimum. The coefficient AL also determines this point. However, this is quite 
difficult to calculate. The approximated value of Me as determined by the 


Newtonian method, does not always provide the necessary precision. 
Let us again turn to the example cited in the previous paragraph. 


In each iteration we carry out several computations of the function F along 
the vector gradient. For this we will assume that de = Shine We will give 


parameter s the values 1/4, 1/2, 3/4, 1, 5/4, and 3/2. Of the six values of 


the function we will select F = F in and the a (values) which correspond to 

this function. We will use this coefficient in future calculations of the 

parameters of the cylindrical bodies. The calculations made in this manner 

are shown in Table 4, 

TABLE 4. 

Parameters imations 20 
othe Nos. of the Approxim ue 
Perturbing a> 
Body 


150 
ty 0 | 427 135 37 ito 
15 ny 08 

a 320 | 317 310 | 308 WO.4 
rh 60 ry i? 6 150 
hy 100 ] 107 Cn ke 
a 570 | 76 | a 590 150 
fy 30 | 137 Wo} 9 200 
f 2 216 | 207 1060 

ee oe 6 


eI 


A — | 0,02102 | 0,024¢ | 0.09137 | 0.0810 | 0.07958 


0.1763 


10728 


Commas indicate decimal points. 


Figure 5 shows graphs of the function change along the vector gradient. 
In the first iterations the minimum of the function corresponds to the value 


Ay = Ary fs = 1). However, right after the second approximation the function F 


acquires its minimum value at a = Shine where at first the parameter s = 3/4, 


then s = 1/2, and s = 1/4. It is evident from curve 7 that the calculated 
value of the function F at s = 1 will be larger than in the previous approxi- 
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mation. Solving the problem with the fixed value s = 1, we also obtained 
sudden jumps in the changes in function F (Figure 5, b). 


. F108 


Fits) 198 


Figure 5. Change in Function F Along the Vector 
Gradient (Various Iterations). 


One must also consider one other peculiarity. After the sharp jumps in the 
value of the function the next iteration, as a rule, is accompanied by a large 
value for the coefficient s. Thus it is evident from Table 4 that, after the 
2lst approximation, a malfunction occurred in the machine (the calculations 
were not being re-checked), the smooth decrease in the function was disturbed, 
and the 22nd approximation became worse than the previous one. During cal- 
culation of the 23rd approximation the value of the coefficient was s = 3/2. /26 
This means that the vector had raised considerably. The change in the function 
along the vector is illustrated in Figure 5c. In following iterations, when 
the value of the function is small, its change along the vector gradient becomes 
less sharp (Figure 5d). 


The example cited indicates that it is necessary to determine the coeffi- 
cient s during the calculation. Henceforth, all computer programs have been 
set up with allowance made for the calculated value of s. 


6. Increasing Convergence of the Fast Descent Method While Solving the Inverse 
Problems of Gravitational Survey 


The nature of the function change along the vector gradient was established 
in the previous section. In this vector, the function F which must be minimized 
is dependent upon one parameter, which is the magnitude of s. The problem lies 
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in finding the value of parameter s at which the function F along the vector 

gradient acquires its maximum value. It is possible to come to a conclusion 

about the nature of this function by using the fast descent geometric inter- 

pretation method. Let us assume that the function F changes according to the 
parabolic law 


ePceas? + bs +0. 
ef eas’ + bs + ¢ (1.25) 


This approximation agrees completely with the computed results shown in 
Figure 5. 


We will find the s at which F (s) = Fil For this purpose it is sufficient 


rn 


to find the root of the equation ae O. By differentiating (1.25), we find 


os 
2as + b = 0 


_b_ 
2a ° 


From this 


In order to compute the unknown parameter s, it is necessary to find co- 
efficients a and b in the equality (1.25). 


We know the value of the function F (0) = Foe At s = 0 the function 


acquires the value of the preceding computation. In order to find coefficients 
a and b, we compute the function F (s) at the two values s = s, ands =S,,. 


1 2 
Let F (s)) = F, and let F (s5) = F Inserting these values into (1.25), we 
obtain 


1 2° 


F =c, 


F, = asi + bs, + Fy 
PF, = as; + 65, + Fy.4 


From the last two equations we find a and b, and after this it is easy to 
determine parameter s. 


5 2a — FH FF) ae! 
7 i 2 [s, (Fy — Fy) — 4 (F, — A) (1.26) 


Thus the following calculation method is recommended. We will designate 
the previously determined value of the function by Fo: We compute function 


F (s) at s = s, and s =s,. We obtain F (s)) = Fy and F (s.) =F We compute 


1 2 2° 
coefficient s according to formula (1.26), and then compute the value of 
function F (s). 


There is a question about the choice of the values S$} and S5- Let 
S$, < S$). The best result in searching for the minimum of the function F (s) 
should be expected when S$) < Sain < 52° In choosing values Sy and $5, it is 
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recommended that one use the value s = Spr? as calculated in the preceding 


iteration, and accept s, = 1/2s__, and S5 = 3/2 of Sor’ In calculations of the 


1 pr 


first iteration, one can assume that $s, = 1, and that S5 SQ, 


7. Another Algorithm For the Minimization of the Function of Many Variables 


We have established that, when solving the inverse problem, the minimization 
is based on the function 


' F=F(p,, Dy +++» PN)- (1.27) 


The algorithm described here is based on the fast descent gradient method 
{[10, 22]. Let us establish some initial value om, oe beats po?) by which /28 
a certain point is fixed in N-dimensional space. Now we choose one direction 
in such a way that it coincides with the vector gradient, but such that its 
direction is the opposite of the vector gradient. If a function increases as 
much as possible along the vector gradient, then in the other direction it will 
decrease, The chosen line is described by the directing cosines 


; —F,, 
| COS @) = —7—= mer. 
Function (1.27) along the vector gradient can be written like the function 


of one variable 2. For this purpose it is necessary to introduce into (1.27) 
the parameters 


Bi =P)” + beosa, Us 1, 2, rie eee (1.28) 
Thus along the selected axis F =F (2). 


Now let us pose the problem of finding the value Z = Z7* for which F (2*) = 
Min. For this purpose it is necessary to solve the equation 


F(j=9() =0. (1.29) 
The new designation of the function y (2) has been introduced for the sake 
of convenience in future work. It is necessary to find the root of the trans- 


cendental equation (1.29). 


We use the reduction method of solving transcendental equations to solve 
differential equations [58]. 


Thus, the equation is given 


y (2) = 0 (1.30) 
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Let us examine the function 
a ised (1.31) 


The value Z = 7* which converts this function to 0, is the root of equation 
(1.29. If function (1.31) has the reverse form 


f= LO, (1.32) 


then the task of finding the root of equation (1.29) is reduced to computing 
function (1.32) at t = 0 because 


The derivative of function (1.32) is like the reverse (1.31) 


- db 1 
ao" OS”) 
Thus we have the differential equation of function (1.32). Having selected 
the arbitrary value Z = Le and having inserted it into (1.31), we obtain the 
initial conditions for solving differential equation (1.34) at t = ty = (20). 
t=, 


oO 


Only one value of function (1.32) 2 = L (0) is of interest to us. 
The integration interval is determined Dy 
At = lino $ini = O— fo = — 9 (h)- 


It is completely natural to select the initial value Ly = 0, and then 
At = —9(0). 


If we use the Runge-Kutta method to calculate Z*, having accepted a com- 
putation step equal to At, then it is necessary to compute 


0 0) =9(0) 70) 
er ea ae 
\2 Pe Lee 


mig: + 2h, + 2p + hy. 


h= 


Then 


Now new values for the unknown magnitudes are determined according to 
(1.28), by 


Pe pO ae co (1.35) 


If, at the computed values of the parameters, the function is sufficiently 
small, then the computation is complete. The formulas (1.35) give the final 
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results of the solution. If not, the iteration cycle is repeated. This point, 
the coordinates of which are calculated according to formulas (1.35), is accepted 
as the initial point. 


8. Program For Solving the Inverse Problems by Using the Minimization Method 


(Fast Descent) 


In its general form, the entire program is written in algorithmic language 
(ALGOL-60). It consists of individual blocks and constructions. The magnitudes 
encountered in the program are described at the start of the program. We will 
subsequently discuss these magnitudes in somewhat more detail. 


In order to solve the inverse problem, numerical information has to be fed 
into the computer. This information consists of several groups. The first 
group is the information which describes the observed gravitational field and 
contains the coordinates of the points and the value of the field at given 
points XT [l:n], YT [1:n], and GNABL [1:n]. The second group contains the 
parameter values of the bodies which approximate the geological model: PP1 
[l1:m], PP2 [l:m], ... PPT [1l:m] and Pl [l:m], P2 [1:m], ... PK [1:m]. In 
massifs PPl, PP2, ... PPT [l:m] the-parameters to be determined (the variable 
parameters) are unified. They belong to functions (1.2) and (1.3). 


Massifs Pl, P2, ... PK [l1:m] determine the geological model but, when the 
problem is solved, they are confirmed and take on the role of constant parameters. 


Since the dimensions of the first and second group can vary when solving 
the various problems, it is necessary to put into the computer memory the /30 
values of the magnitudes n, i.e., the number of points which are used in 
solving problems m, i.e., the number of perturbing bodies through the use of 
which the geological model is approximated. With regards to the magnitudes 
(T + K), i.e., the general number of parameters in an elementary geological 
body, this is constant for every actual problem. For instance, if a geological 
body is a made up of cylindrical bodies, then T + K = 4, T = 3 (these are 
parameters t, h, d) and parameter K = 1 is sign o. If the geological model is 
described as a group of straight projections finite in space (Figure 1), then 
most often T = 1 (this parameter is d, the position of the projection's 
boundaries). The magnitude K = 5 (0, h, H, Ly» 14); these are all considered 


constant and are not subject to change during the minimization process (1.2). 


In addition, values Fein and Acin must be placed in the machine's memory. 
These magnitudes determine the end of the computation cycle according to 
criteria (1.13) and (1.14). 


Value Feo, can be computed, but magnitude Agia 
must be determined during the process of sampling the actual practical problem, 
or it must be established from research with models. Most often this magnitude 
varies within the limits of 0.05 - 0.02. Thus, the massif of the initial value 
is determined within the operational memory device in order to solve the problem. 
We will designate this M,- The dimensions of this massif depend on magnitudes n 
and m. 


is assigned. This magnitude 
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Data fed into the computer is printed. This is necessary because one might 
want to check the initial data. After all, even an error in the numerical data 
fed in which at first appears insignificant can distort the intent of the in- 
terpreter, and the problem will not be solved. 


Now we again turn to a description of the variable magnitudes in the 
program. The variable indices are designated by the symbols i and j: i changes 
from 0 to n, and j changes from 0 to m. Magnitudes SKL, TKL, P35, and P7 are 
determined within the program itself. The first two magnitudes are needed to 
transmit directions to switches KLM2: = AO, Al, A2, A3; switch KL37: = T6, T14, 
D. Near the switch designation, marks are indicated by means of which several 
instructions in the program are designated. The magnitude SKL, during the 
process of solving the problem, can acquire a value from 1 to 4. Depending on 
the value of this magnitude, switch KLM2 sends directions to one of the 
operators with mark AO, Al, A2, or A3. Magnitude TKL can acquire values from 1 
to 3 during the calculation process. The entire numbers P35 and P7 are also 
used to transmit directions. 


The massifs of the operative numbers are then described. First come the 
coordinates of the points and the values of the field at the given points, and 
then come the parameters of the geological model. In order to record the in- 
termediate data, two groups of massifs, M, and Mz> are distinguished. Within 


group M, massifs PP] and 2, PP2 and 2, ... PPTM2 [1:m] are distinguished in 


order to record the results of the calculations in the regular iteration of the 
geological model's parameters. In each cycle of the calculation process, after 


changes have been made in the geological model, the gravitational field Vine 


(x; y;) is computed at n points. In the program this massif is designated as 


GTM2 [l:n]. The calculated values of the geological parameters and the cal- 
culated gravitational field (massifs PPl and 3, PP2M3, ... PPTM2 [1:n] and 
GTM3 [l1:n]) are re-written into massif M3. Symbols FPR1, FPR2, ... FPRT [1:m] 
designate the massif of the values derived from function (1.2) according to the 
desired parameters. 


Then follows a description of the operative number. S has already been 
indicated, Fein and bein are introduced at the same time as initial data. They 
are described by the identifiers FKON and DELT. We will discuss the value of 
the other numbers later. 


Now let us go to a description of the program itself. The values corre- 


sponding to 1.0 and 2.0 assume magnitudes $y and S5- The whole number P7 


receives the value 0. Then follows the component operator, described by the 


Sign M. The values of the variable parameters in massif M, are rewritten within 


the cycle. Variable P35, SKL, and TKL acquire the value 1, and directions are 
given to the operator assigned the symbol A. 
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At the very end of the program are located four groups of operators, next 
to which are the symbols A, B, C, and K. They are closely interrelated and 
are essentially specialized sub-programs. 


The symbol A designates the block in which the values are computing 


i theo* 4) = ps Vitheol Xi Yi). i 


any oe p> [Vopse (xis yi) — “theo (xis yd)? 


The block operates four times in one computer cycle. Switch KLM2 (switch KLM2: 
= AO, Al, A2, A3;) is located here. The function F is written in tur into 
four different memory cells. 


The next block (symbol C) is activated only if the number P35 is not equal 
to 0. Here the derivatives of the function F are computed according to the 
desired parameters. In each cycle these magnitudes are calculated only once. 
This means that every time, when switching to A (and this means to block B as 
well which follows it), the values of number P35 must always be noted. 


Using operator K, an output is made from the unique sub-program into the 
necessary position in the basic program. 


It is worth noting that operators A and C cannot be actually described in 


the general program. In each case its function will be Ms theG (x, yj); and 


this means that they will be various values of the derivatives of function F. 
When examining the actual problems, it is necessary to describe these blocks. 


One moves from the sub-program with the symbol A - B - C - K to operator 
T6 at (TKL = 1). The computed value of function F, the massif Ves (x; 5 y;,); 


and the values of the variable parameters are rewritten into massif Mee 


Directions are then given to operator T9. The decision is then made to print 
the direct problem. Thus, there is the possibility of comparing Vee (x, y;) 
with Vanes (X55 y;) computed for the model of the first approximation. Here 
the entire computation process can be stopped, and a corrective factor can be 
introduced into the model of the first approximation of geological structure if 
the interpreter decides that the correction was not constructed satisfactorily 
earlier. By means of operator P7: = 1, access to operator T9 is closed in the 
future. 


In block T10 the coefficient AEN is calculated. Then (operator T11l) one 


determines the values of the desired parameters for s = Sy and s = So» while 


each time reference is made to the sub-program (blocks A - K). Before moving 
to block A, parameters P35 and TKL are given values corresponding to 0 and 2. 
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This is necessary in order to avoid computing the derivatives (block C), and to 
return to operator Tl4. Here F = Fin is searched for, and the values of the 
parameters in the geological model, as well as the massif Viens (x; 5 Y;) which 
correspond to them, are rewritten into massif Mz. 
The coefficient s is then computed according to formula (1.26), and the 
values of the geological parameters are again computed (during this process 
instructions are again transmitted to sub-programs A - K). Symbol D designates 
the operators which minimalize the calculated function F. Some relationships 


are then checked. If Fein (which has identifier FM3) is equal to Fy then this 


is evidence of the fact that the function has not decreased in the given cycle. 
We have reached the possible F as for concrete geological data. The computation 
should be finished. 


The end of the computation occurs when eres is less than or equal to Pein 
or when the character of the function's monotonic progression is sufficiently 
small. Meanwhile the results of the minimization are printed out: the values 


of parameters, the collected function, and the fumction F itself. (symbol D1). 


If none of the described criteria have been fulfilled, then it prepares to 
move on to a new iteration and directions are transmitted to blocks N. 


Minimization Program for Functionals 


.integer n, m, I, J].-P7, P35, SKL, TKL, T; 
-real array XT, YT, GNABL (I:n], 
; PPt, PP2, ... PPT [!: ml, 


—— 
Begin 


: Pl, P2, ... PK [1:m), 
t Tha os REID, ae . PPTM2 {1:ml, 
_ 
aha PPIM3, PP2M3, ... PPTM3 [1:m], 
GTM3 {1:n 


FPRI, FPR?, . -FPRT [L:m];_ - 
4 real FKON, DELT, FM2, aye Fl, F2, FF, 
FM3, LAM, si, 82, S$ 
switch KLM2:=AO, Al, A2, Maa; j 
éwitch ag me ce D; mw 


Input — 3 


=" array xT, af. ‘GNAB fakes o 

oP TU mi), = 

5 Ph, a ” PK my} en 4 
Printout Begin real. ppt bP, XT, YT. ONABL’ (ised, . 

PPT [\i:m], -. pe 

BL ‘PK he ml] end 


Begin . tual 0; $2: 30; hile end 
N: Begin fotf=T step.1 unti 
Begin gat tl iPr PPIMOLj|:= PALI feat 


jJ en 
: P35: Mr =TKL:=1; go to A end 
ae "Begin PPIMGLL wre IMI: PP2MSJI: =PP2M2Ijh 
“=Begin = = 
pPIMa Maijl-=PPTMaL] en 


’ for eh el ae 1 until n do 
GTMMT]:=GTMiI; 
a if P7 = 0 then go to T9 else go to T10; 
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gin: P7:= 1}; 


ps Beg ntout real GTM2 [I: nl end 
Te ti out sl EOimrny } 


ae cere, faa 0.0 
‘ 
* 


ee ont step 
Begin: K Wiel ile FPR@&1t 2... + FPRTI ¢ 2 
dane =sum+KW 
vee =EM2/sum_ ae 
| ur tor Sw: pean 
Par we 1 until m d 
P1M3 UE = PP] ul _ sw x LAM X FPRI[j};~ 


oer a Se eo. 


PPT Nat = PT r= SW <x LAM X FPRT[j] end 
P35:=0; TKL:=2; go to A; 
Tt of FM3 > Hoc art we 
* Begin - FM3:- FM2; | 


. for fed step | until m do * | 
2 Beeln PLM ifm PPM Ul x 


vibra | 


t, 


“PPTMS = PPTM2 ted’ “3 
a ghar trad te oe n do 4 
wee Soke ei 

} i, 1p aX (FL — aay FE x (F2-FO)) 


for 1:1 ster 1 until n do 
‘  GTMAlil: = GTM}; 
» for j:=1 step 1 until m do 7 
ae Pon PP1M3{]]:=PP1M21)); se 
PPM: Eerie | 
‘M3 = FO) V (FMS « < FKON) V (abs (((FM3 — FO)/FM3) 


2 T)) then go to Di; 
Begin see leg eeee pigs directions to printout array GTM3 fi: nj, /34 


for j:==1 step 1 until m do aac: a . 
Begin PPIIj}:= PPIMAIil; ie: 
cet PPTI|I:= PPTMAL end 7 


et Se =0.5 x SW; S2:=1.5 X SW; go to N; 
Di:Printout eLFMS: fi GTM3 [1:nl, : 


1 .-. PPTM3 11: Stop! ha, 
mij ll: im), description V1; a 


| + As Begin real lie 
wa: for 1:=1 mf T entil n do ae F 
f-a% Begin GTM2 [1}:=0.0 
eee for j:=4.step | until m do 
. Begin ay Ly == Vi (PP1M2, ... PPTM2, PI, oy YT): 
1); = ore sil + DGTIj} end end *' 

B: _.. Begin ra F, Olea ae 

he: = 0.0; ; 


he for |:=1 step 1 untl! n do 
“: Begin RAZN:=GNABL{i] — GTM2 [i]; 
RAZN2:=RAZN X RAZN; 
F:=F + RAZN2 end 
es FM2:= a to KLM2 (SKLI; 
_ AO: F0:==F; go to LI; 
. Al: Fl:=F; go to L1; 
es AQ: FoiaeF: go to LI; 
ae _, Ad: FF:anP; 
“ye. LE: SKLe SKL + J; 
a. 5° If P36.== Q then go to K end 
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C: “Begin real ec Lf oT t =O FPRT:, 
for j:=1 step :J- - es — 
Begin FPRIER a 
for i:=1 step. 1 tail do 
Begin f1li]:=f - 


% FPRTI/jI: —FPRit + {T[i] end end end 
: Ke o ee KL37|TKL} end end 


OO heer eae es 
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CHAPTER II 


ESTIMATING THE SIZE OF PERTURBING BODIES' EXCESS 
MASSES AND CENTERS OF GRAVITY 


1. Solution of the Inverse Problem of an Anomaly in Gravitational Force for 
a Group of Spherical Perturbing Bodies 


Let an anomaly in the force of gravity be given. Let us assume that the 
perturbing masses are isolated bodies. Some of them can be located relatively 
close to one another, thus creating a general, complicated picture of a gravity 
field upon the surface. In order to estimate approximately the location and 
dimensions of the perturbing masses, we will identify their form as spheres. 
The position of each sphere can be characterized by these parameters: X52 Yo 


are the coordinates of the epicenter, h is the depth to the center, and M is 
their excess mass. 


In an observed field it is almost always possible to establish a number 
of perturbing objects. We will designate by M the number of spheres by which 
the perturbing bodies are identified. Now the observed gravitational field 
on plane xOy can be approximated using the following formula: 


~~ 


Aptech Ya (2.1) 


ans 


Is is most convenient to express linear values in kilometers. We will 
consider the excess mass in units of 109m. For instance, if the excess mass 
is 12° 1010p, then the formula should be written as M = 120 units. When 
selecting units of measurement in this way, the coefficient k in formula 
(2.1) must be accepted as equal to 6.67, and the gravitational force anomaly 
will then be expressed in milligals. 


Each sphere is designated by four parameters. Four M parameters go to 
make up expression (2.1). 


In the observed field we will fix N points with coordinates (x; , y,)- 


We will include in the number of these points the most descriptive points of 
the gravitational anomaly. We will establish the function 


=—— ws -— 
oe ee eee - 


Fa 2 (AB obsel%e> Ys) — Ag eo %i» Had" (2.2) 


The function 8 heo (s, y) is determined by formula (2.1). Relation- 


ship (2.2) contains four M parameters (xQ5 Vag? h., M.), i ig 22 Sia eat 


J 
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Of the total sum of the proximate values of these parameters, we will /36 
choose those which can convert function (2.2) into a minimum. We will mini- 
mize function (2.2) by using the fast descent method. All computations are 
reduced to a sequential determination of parameters according to the formulas 


ft? = a (Pale 
| it = yy =, oe (Fa, )e» “ 7 ak 4 
ee ee 4 
SMP = My Fe 


(2.3) 


where k is the iteration number. 


The computation model was described in sufficient detail in the previous 
chapter. Now we may cite only values for derivatives of functions (2.2) 
according to the pest ted parameters: 

pe 


bate, igh. 
Figg, = 4: BRM A y 
Oj = 6b atlas — a + es bel" + 


_. 6d, SYS = tl) s 
Foy = rp) (ler — x — i = sh fs 


2.4 
ym — ag, Sy a el wd, ee 
Ye 7 fed (x; + rg? + (yi — y+ AB” 


6; 
oe) ee ne 
se 1 [ap — ag ty. — oe) hi) i 


where-65 = Ag ee is V3) — Ab greg a? Ya? 

It should be noted that the inverse problem for a group of spherical 
bodies leads us to convex programming. This follows from the theorem of 
V. P. Zidarov [116]. Let the distribution of the anomalous field of gravi- 
tational for ce be known on plane x0y (V = 0), and let it be established 
that the anomalous geological objects can be approximated by N spherical 
bodies, each of which has excess mass ML: 


Let us assume that the centers of these masses are concentrated at 
points Q (Ey. "e oy): We examine the ete ered 


where A” me ' is the value of the gravitational force from the spherical 
‘ ge RIE) /37 
keel ed | —_—_— 


perturbing bodies at point P (x, y, 0) of plane S, and t, a= Vi =, +e Ns)" +e 
is the distance between points P and Q.. 
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It is proven that there exists only one position of points K_ at which 
U has its minimun. 4 


2. An Example of Solving the Inverse Problem 


Figure 6 shows a gravitational force anomaly. Here four perturbing 
bodies are quite clearly distinguished. The position and size of the excess 
mass can be characterized by 16 parameters. Let us establish the beginning 
of the coordinates at some certain point in the field. We will superimpose 
coordinate axes 0x and Oy on the plane where the anomaly has been determined. 


Figure 6. A Gravitational Force Field Caused 
by Several Isometric Masses. 


We will confirm N points which would best describe the observed field. 
We include in the number of these points the field's extreme points, and some 
points in the area where the field changes monotonically. These points are 
noted by circles. Numbers have been placed next. to them, indicating the 
observed value of the anomaly and the coordinates of the point. We will 
record the coordinates of these points and their corresponding gravitational 
force anomalies (Table 5). In all 31 points (n = 31) have been fixed. 
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TABLE 5. 


A, SCT OE NG Oy Aa 2 ee 
ve 4 ich 5 MMe tie ~ FUN Ay “3 
Coor- fame ™ 7) Coor- a ie 
A dinates : g 5 # i : dinates 
26 Gs ts oyltor= : 
JE oa 14 | Ba 
oof kaa ae ae WZ “O ' 
1 0 0. 7,73 F 
2 1,5] 0 68, 
3 | 25]0 2,62 © *4 1.0L 1,6 
4 4,5 | 0 0,48 
5 O |—1,0] 3.19 
6 1,0 J—1,0] 3,08 - 
7. 3.5 J—1,0] 1,19°1 
8 1 1,0-j—2,0) 3,14 
9...) 2,0 |—2,0| 7.5 
10" }-2B.|—3,5) 1.37. 
ll 0 0,32 ; ag 


Commas indicate decimal points. 


In the observed field we will determine the initial values of the per- 
turbing bodies' parameters. We will establish the coordinates of the per- 
turbing bodies' epicenters on the basis of the field's strongest points, 
Then we will examine the anomaly from each point individually. We will con- 
struct a graph of the field, according to its most descriptive profile. We 
will determine the depth to the center of the masses according to the well- 
known relationship h = 1.334 7/2? where X1/2 is the distance from the point 


of the field's maximum to the point where Ag = 1/2 Ag ax’ It is possible to 


determine the size of the excess mass according to the relationship 
M = 150 h@ Ag ax’ Here the anomaly is expressed in milligals, the depth is 


expressed in kilometers, and the excess mass is expressed in millions of tons. 
In our example the initial approximations were deliberately made worse than 
those obtained through computation. We will reduce the geological premises 
into a single table (Table 6). We will not compute the value of the function 


F at which the computation should be ended. We will make use of formula (1.13). 


We will assume the observation error to be equal to 0.25 milligals. In this 


case ) 
F = 2 ¢ 31° 0.0625 = 3.88 mgl. 


In Table 7 we have shown the change 
of the minimized function from itera- 
tion to iteration. The computations 
continued until the seventh approxi- 
mation where F_ = 2.99 mgl“, This 
corresponds to an average deviation 
of the calculated field from the ob- 
served field of 0.23 mgl. 


Commas indicate decimal points. 
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TABLE 7. 


ar ; cme, Sma ¥ ay 5 a np ey a “4 
| Number of . 
: Approxi- 
12,7 


F 225. 172 | 15,6 


9,33 | ss 4.55 | 299 | 


Commas indicate decimal points. 


In Table 8 the results of the calculations are shown. Since a theoretical 
example was examined, it is possible to compare the observed parameters with /40 


those which were assigned when computing the anomaly. The latter are given 
in Table 9. 


___.__. TABLE 8. TABLE 9, 
: ee Seenee eee F 


Number of the 
Perturbing Body 


a 


| 0,58 | $ 
179-10" | n|ter. OF 0.135.10°| | 
Commas indicate decimal points. Commas indicate decimal points. 


In Table 10, the gravitational force field selected at the 31st point 
is compared with the observed values. Here are also given the values of an 
anomalous field from bodies of the first approximation. 


TABLE 10. 
2 } Coor- Value of Ag= 
wo i eG 
fe 32 dinate | 3 &| 
a a we a > D ore 
‘2% é A 
, | c | ee 
1Jo fo each 1 4,8311.0 
271.5 10 3,86 . ‘2,0°T 2,0 : 
3/25 | 0 2.98 3,0 | 0,5 
4}45 10 0,47 ) 4.0 | 4,0 
5 | 0 —1,0 33 ) >? 2,0 |—1,0 
6] 1,0 |—1,0] 2,70 3,03. | 3, 2. | 2,5 |—1.0 
7)35 |—1,0} 0,89 1,15 + 3 | 2,0-] 0 
8] 10 |—2,0] 2,66 3,12 3.14 7 —1,07 0, 
9} 2.0 |—2.0| 7,33 7,58 7,57 [25 | 2,5 |—2,5} 
10} 2,0 |~3,5} 1,07 1,35 1,37 | 2] 2.0 |—2,5 
1] [—3,0] 0 0,24 0,31 0,32 | 27°F 3,0 3h. ‘ 
12 |—1,0/—1,0{ 1,37 1,66 1,67 2.0 | 3,0: 
13 |—0,5 |—2,0| 0,81 1,02 rae .1 0 1,0 R , 
14:1—1,0] 10°] 1,38 1,67 1.6 . 1,0] 0 G1.| 3,79 | 3,82. 
1510 25 0,80 oat 1,03 5 31 | 4,5 |—4,0| 0,18 0,24 Dee 
16/05 {1.5 | 191 | 29%) 2:35 


Commas indicate decimal points. 
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3. A Practical Example in Evaluating the Depth of Location and Mass of 


a Group of Perturbing Bodies 


In one of the areas, gravitational and magnetic survey sketches were 
made. After excluding the regional component, a sufficiently clear gravi- 
tational maximum (Figure 7 a) was distinguished. 


eae a : a 


Figure 7. The Gravitational Force Field and Magnetic Field 
in One Area of the Investigation. 


In this same area three local objects (Figure 7 b) were distinguished in the 
magnetic field. One can make various conjectures about which perturbing 
objects caused these anomalies. It is possible that the bodies which cause 
Magnetic and gravitational anomalies cannot be compared to one another. 

In a given region rocks can exist which create both magnetic and gravitational 
anomalies. There can be any number of intermediate variations in which dense 
rocks are distinguished from magnetic rocks, or there can be rocks which are 
characterized by excess density and by higher values of magnetic properties. 


The researcher who is thoroughly familiar with the geological structure /41 
of the region in which he is operating and with the physical properties of 
the rocks can choose one hypothesis or another. 


The following problem was posed to us. On the basis of the geological 
premises in a given region, it is possible to propose the existence of four 
perturbing objects with excess masses. It is necessary to find four centers 
of mass and to determine the excess mass of each body. 


In order to solve this question we approximate the perturbing masses in 
the very first approximation by using four spheres. We will fix 22 points 
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(n = 22), choosing them in such a way that they describe the field being in- 
terpreted as completely as possible. Table 11 shows a listing of these points 
where the coordinates and values are given for the gravitational force anomaly. 


We calculate Fen according to formula (1.13). We set S6Ag = 0.2 mgl, and 


a 2 
then Fein = 1.76 mg1 3 


TABLE 11. 


Coor- SPs C aoe, | 
_ dinates [ A i - Value of & Bae | 


‘Ob- 
Serv- 
d 


‘oe 
» 
> 


1. 142 1 

a 46 ‘142 23 | 23 

3 &6 | 52 2, 40 

4 43 | 40 2, 3,0 
8 B2 | 53 21 | 20 

7 44 145 0, 

8 4,7 

9 1,7 

0 23 


MRONO Vee NW 
oS 
i) 


ooro= 


~~ Commas indicate decimal points. 


Table 12 shows the parameters of the four spheres by which we approxi- 
mated the perturbing bodies in the first approximation. The results of solv- /42 
ing the problem are shown in the same table: in Table 12 are the parameters 
of the selected bodies, and in Table 11 the observed and selected values of 
the anomaly have been compared. As is evident, the selection was made with 
sufficient precision. The maximum deviation is observed at 3 points and is 
0.4 MGL. 


TABLE 12. 


|e of Perturbing Body 
: Parameter ©) |}——#_{———_—_—_ 
1 fie be a |} 


Initial Approximation 


po x “ea | 075 | 085 | 15 
' y 0 0,55 | —0,15 0,4 

h 0,3 0,35 0,25 0,3 
tn Units-of M+10°m | 0,19 0.20 0,18 0.17 

*y «Solution 
x 0,38- 0,79 0,93 144.8 
y —0,02 0,62. | —0.20 037° = 
h irs 055 0,57 0.55 

In Units of M+ 10°m| 0,078 | 0,124 |. O18 0.139 


Commas indicate decimal points. 
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4. Solving the Inverse Problem for Anomaly vies for a Group of Cylindrical 
Perturbing Bodies 


The solution to this problem has essentially been given in the first 
chapter. Here we will simply explain what concerns the immediate problem and 
we will make some remarks on methodology in selection of the initial approxi- 
mation. 


The function which must be minimized is written by relationship (1.23). 
Formulas (1.24) are used to write the recurrent relationship to determine the 
desired parameters {t,, hj. ai}, j = 1, 2, ..., m, where m is the number of 


cylindrical bodies. 


Now let us examine the question of selecting component values for the 


first approximation: oe bi a 


% 


» sign oF) (all information about the 


geological structure must be taken into account here, or hypotheses about the 
distribution of the perturbing masses must be introduced). The parametric 
values of the cylindrical masses for the first approximation can be computed. 


At the very start it is convenient to determine components gto). gl) 
can be accepted as equal to the values of the abscissas of the point where 
Me = 0, or the points relative to which the function is symmetrical in a 
narrow region. Each local anomaly can be examined individually. By fixing 
the beginning of the coordinates at points x, = o each time, we computed 


the initial approximations of depth according to the relationship he = x V3 
(x, is the abscissa of the extreme point). The initial values of the para- 


meters can be computed according to the known relationships for the mass of 
a cylinder 


Oya a 
ay = a aw (V,;) mary 


fe ®t 


or gk 
: Ap. 0. ol 1h} (V,.) max: , 
Soe. ES a 
If the magnitudes h are expressed in meters, and if er is expressed in 
Eoetvoes then x. will have the dimensionality m/m. Now it is easy to cal- 
culate the parameter 160) =vA.,. 
j i/p 

In conclusion, let us examine a practical example. In one of the ultra- 

basic massifs, an anomaly was obtained for the horizontal gradient of the 


force of gravity (Figure 8). By their nature the perturbing masses can be 
identified with two-dimensional bodies. A detailed examination makes it 
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possible to distinguish five perturbing masses. One body has a very deep 
location, and its anomaly has determined the general nature of the field. 

The four other bodies are small masses, and their location is not deep. They 
cause anomalies in the background of the basic field. We will identify each 
perturbing body as a cylinder. Since the position and dimensions of the 
cylinder are characterized by three parameters, 15 magnitudes form the basis 
of the determination. By taking into account the nature of the fields, we 
will fix 17 points equidistant from one another. Table 13 shows the numerical 


values of the abscissas of the points and the anomalies which correspond to 
then. 


nee PETES eR 
£ ae | 
ae 


{a Eoetvoes be oe Si UME o | 
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: : ey ee 
Figure 8. Anomaly V,, in One Area of the Investigation. 


TABLE 13. 


\No. of | ° 
"the point/*» ™ x, m No Selecte 


SF 


Ox nOnpotes — 


On the basis of the observed field, we calculate the initial approxi- 
Mations of the unknown parameters. These are shown in Table 14, We will 
determine Fein according to formula (1.13). We will accept 6v = 5 Eoetvoes, 


and then Fein = 850 Foceioss Tables 13 and 14 show the results of the 
solution. The initial value of the minimized function F is 4399 Foetioes /44 
and the final value of F er is 622 Eoetvoes*. This corresponds to an average 
deviation of the observed and selected fields of within 4.2 Eoetvoes. As is 
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evident, the selection has been made with sufficient precision. In Table 13 


the observed and selected anomalies have been compared. The maximum devia- 
tion between them is 13 Eoetvoes. 


TABLE 14. 
Bes : i om SS - 
Parameter of the Perturbin Bod - 


we, age Initial Approximation bs 


35 , 
620 
ee 
avi this fie he ees 
“MT OY 23 23 | 3.9 
162 - 75 85 12.2 
196 425 526 633 
6,9-10° 1.7-108 1,7- 108 0,05- 103 


A he eee 
‘. ——- sends ses ee —-— ee, S28 Pj Bena 


Commas indicate decimal points. 


CHAPTER III 


DETERMINING THE CONTOURS OF TWO-DIMENSIONAL GEOLOGICAL 
BODIES ON THE BASIS OF GRAVITATIONAL ANOMALIES 


In this chapter an examination is made of the questions of solving inverse 
problems of gravitational force and horizontal gradient anomalies. The unknown 
controu of the perturbing body is approximated by segmented straight lines with 
right angles. 


The methods of solving the problem are illustrated by theoretical and 
practical examples. 


1. Solving the Inverse Problem For Two-Dimensional Perturbing Bodies Bounded 
by Segmented Straight-Line Contours 


On the basis of the observed gravitational anomaly, it is necessary to 
determine the contours within which the excess masses are located. We approxi- 
mate the unknown controu by using a segmented straight line. In this case it 
is possible to consider the gravitational anomaly as the sum of the effects 
from oblique or straight steps. We will accept that within each step the excess 
density is a constant and known value. By changing the individual density from 
step to step, it is possible to describe the configuration of a geological body 
with a variable density. 


In order to simplify the solution somewhat, we will give the segmented 
Straight-line contours right angles. If there are m steps, then the anomaly 
can be approximated in the following manner: 


For the gravitational force anomaly 


m = 
Ag (x) ke Voo;|atH, — hy + 2H; arctg ai i _ (3.1) 


=l\ i 


7eng sare 


—s 


— 2h; arctg —— 
! 


dy’ ae Hi +(x — dj 
Pe eat 
Gee a eee? Lae 


and for the horizontal gradient anomaly of the force of gravity 


Le ae ae wa iE 
S HY + («— dj)? 3.1a) 
=k Yon te 
ce wa "RR EG—ap 


Each projection is characterized by four parameters: he and a which are 


the depths to the upper and lower boundaries of the projections, a which is 
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the abscissa which determines the location of the vertical boundary relative to 
the beginning of the coordinates, and a which is the excess density of the 


masses. The value of these parameters was shown in Figure 1. In this case a 
two-dimensional problem is being examined where L > - © and L, + + ©, 
We will fix in the observed field n points with abscissa X;- We include 


those points which are most descriptive of the gravitational anomaly. We will 
construct a function like (1.2) which will now be expressed in the following 
manner: 


For the gravitational force anomaly 


Pe x [A gobse(*:) = Agtheo(XOl (3.2) 
For the anomaly of the horizontal gradient of the force of gravity 
(3.2a) 


i : 
jF= p> Vezopse(%) — Vertheo (ol : 
eet ve 


The values of Lee (x) and V (x) are determined according to formula 


(3.1) and (3.1a). 


xz theo 


Thus the function F will contain four m parameters: o,, Nes are and qd 


7 
(j = 1, 2, ..., m). We will assign the values of parameters h., a and o.. 
di; remains the variable parameter. Of all the possible values for these para- 


meters, we Will choose those which may change function (3.2) into a minimum. 
We will carry out the calculations according to the formulas 


dt? = dy — Ay (F'di)e 
where k is the number of the approximation (iteration). In these expressions 
mG acquires the values: 
For the gravitational force anomaly  _ 
= J20kg; NG ln eo eae, 
Fay ARO OE an aie (3.3) 


8 = Agobseti  Agiheol(tie 


And for the anomaly of the horizontal gradient of gravitational force 


cs xi — 4) (5 =f) 
Fay = — 40 918: 5 ela 


St Sa EP aan 2 3.3 
HEL xr — 4p MAE xi — a (5-32) 


8, = Vaz obse (Xt) — Viz theol*0)- - 
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In these formulas k is the gravitational constant. It is most convenient 
to express linear units in kilometers. In order to express the gravitational 
force anomaly in milligals, it is necessary to accept the coefficient k as 
equal to 6.67 in formulas (3.1) and (3.3). When calculating the anomaly of the 
horizontal gradient of gravitational force, it is necessary in the formulas {47 
(3.la) and (3.3a) to accept k = 66.7, whatever the dimensionality of the linear 
values. In this case the anomaly is expressed in eoetvoes. 


2. Method of Selecting the Contours of a Perturbing Body 


We will assume that a gravitational anomaly is given and that it is 
established that the perturbing body is two-dimensional. Let us examine how 
the interpretation is done by the selection method, using the measuring grids 
of Barton. All half-space z > 0 is divided up by straight lines z = n into a 


number of elementary strips. In each of these strips, right-angle areas are 
distinguished by straight lines x = d. (t =1, 2, ...). It is assumed that the 


perturbing masses are concentrated within these areas. 


gr oteee do in all 


the strips such that the anomaly calculated from the selected body coincides 
sufficiently well with the observed anomaly. 


The problem consists of finding a series of numbers d.; d 


Let us assume that we know some value z = ho» beginning from which the 
areas with excess density can be found in half-space z > 0. We will fix a 
certain sequence th} with a finite number of digits. 

The magnitudes ae can be determined according to the following relationship: 
as has been accepted in the Barton measuring grid. In this case it is necessary 
to assign (in addition to ho) the values h, and 6. The relationship (3.4) 


determines the thickness of the strip of some geometric progression with the 
denominator 6. In general it is not necessary to compute the sequence th}, but 


rather it can be assigned either completely arbitrarily or by using the geological 
model. 

Thus the finite sequence tht is given which forms m strips. In each j strip, 
the group of numbers td) 5» dj) seus des} are determined. The number of digits 


of each group can vary. These numbers describe the location of right-angle 
regions with excess masses. At an arbitrary point on axis Ox these masses create 
an anomaly which is determined by the equality (3.1) or 3.la). Now it remains 

to fix n points and to establish function (3.2) or (3.2a). 


-, dj. 


The magnitude F is a function of the parameters: F = F (dj; d,; a0 p 
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We will illustrate the solution of the problem by the following example. 
Let a gravitational force (anomaly (Figure 9) be given. We will further assume 
that the possibility of classifying the perturbing bodies as two-dimensional 
has been established. We will choose the system of coordinates in the following 
way. We will direct the axis Ox along the observation profile. We will fix 
the beginning of the coordinates at some point in the profile. We will direct 
the axis Oz vertically downward. 


49 selected aa 
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The Selected Cross Section 


The Initial Approximation 


Figure 9. Example of the Selection of a Geological 
Cross-Section (Two-Dimensional Case) For a Gravi- 
tational Force Anomaly. 


On the basis of general information about the geological structure of the 
region and taking into account the anomalous field, a geological cross-section, 
i.e., the first approximation, has been constructed. The surface layer is a 
homogeneous rock mass throughout its thickness. The thickness of the rock 
begins to vary at a depth of 10 meters. We will make a column graph like a 
Barton measuring grid, while taking into account the constructed geological 


model. The first line z = ho = 0.01 kilometers identifies the boundary between 
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the surface layer and the bedrock. Up to a depth of 1 kilometer, four more 
horizontal strips have been made, each of which is 250 meters thick. The 
thicknesses of the following Sth, 6th, and 7th strips have been accepted as 
equal to 500, 750, and 1250 meters, respectively. This makes it possible to 
describe in considerable detail the configurations of the geological bodies by 


using segmented straight-line contours with vertical boundaries. In this problem 


23 projections have been distinguished. Their location is described by the 
selection of four-dimensional vectors Noss hss Hes dt. The values of these 


parameters are shown in Table 15. 


TABLE 15. 
No. of the body 


Parameters of the 
perturbing body 


Qo +-0,04 |+-0,04 | +-0,04 | +0,04 | —0,04 | —0,04 | —0,04 |—0,04 
h 0,01] 0,25} 05 [* 0,75] 0,01 0,25 0,5 0,75 
H 0,25] 0,5 0,75 | 10 0,25 0,5 0,75 | 1,0 
d 0,75| 9:75; 0,76:4 0,75 1,75 1,75 1,75 1,75 


Parameters of the 
perturbing body 


0a 
h 
H 0,25] 0,5°-42. 
d 40} 40740 | 40 | 40 


4.075 | “Lo £5 2. 
4 


' 


Parameters of the 
perturbing body 


No. of the bod 


a Tra 


Commas indicate decimal points. 


Let us now turn to the anomalous field of gravitational force. On curve 
Ag we will distinguish the points by which it would be possible to fully esta- 
blish an anomaly, according to which the geological cross-section will be 
selected. A list of these points is shown in Table 16. Tables 15 and 16 
present the initial data for solving the problem. We will also establish the 
value F of which it is possible to conclude the computation. We will make use 
of formula (1.13). We will assume that the selection error should be 0.25 MGL, 
and then F,._ = 2.43 NGL’. 

fin 

Figure 9 shows the results of the solution in graphic form. As was said 

earlier, the solution of the direct problem can serve as the criterion for the 


correctness of the model of the first approximation. The gravitational effect 
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from the geological model hypothesis is obtained as an intermediate step in the 
computation. This effect is constantly cited when solving any problem. These 
results can be of considerable help in analyzing the final solution. 


TABLE 16. 


Jokes ets 
Coordination of nt 
the anomaly 


Coordination of , 
the anomaly : 


No. of ‘the poi 


ugh ne E 


Commas indicate decimal points. 


Now let us turn to solving the problem concerning the anomaly Noe Figure 
10 shows a graph of an anomalous field. In the observed field Weg between the 


1.5 kilometer and 3.0 kilometer peaks, a geological object is distinguished, 
the rocks of which possess greater density than the enclosing medium. The 
nature of the anomalous curve in the area of the 2.0 - 2.5 kilometer peaks is 
evidence of the heterogeneity of this object. In addition, in the right part 
of the profile, contact between rocks of various densities is clearly distin- 
guished. The location of the contact is approximately described by the maximum 
of the anomaly V__. 
XZ 
Taking into account data on the geological structure of the region and the 
physical properties of the rocks, and subjecting the observed anomaly to 
stringent analysis, we create a schematic hypothesis of the distribution of the 
geological object. Following this schematic, we assume that the upper covering 
over the entire area is homogeneous to a depth of 250 meters. The lower boundary 
of the rocks, the density of which is 3.4 g/cm, is located at a depth of 1.5 
kilometers. Within this rock mass, rocks are distinguished with a density of 


2.9 g/cm. Their lower limit reaches 850 meters. 


Beginning from the peak at 4.75 kilometers, a second rock massif is dis- 


tinguished with a density of 2.9 g/cm’. On the basis of the nature of the field, 
it is possible to propose that it is relatively homogeneous. The lower boundary 
of the contact, on the basis of geological data, has been accepted as equal to 
3.0 kilometers. The enclosing rocks over the entire expanse are homogeneous, 


and we accept their density as equal to 2.6 g/cm. However, between the two 
massifs described, the enclosing rocks have undergone changes, and their density 


should be accepted as equal to 2.5 g/cm’. Having made all the constructions, 
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we divide the lower half-space (z > 0) into several strips as is done in the 
Barton measuring grid. The first line should be drawn at a depth of 0.25 
kilometers. It will divide the upper covering from the heterogeneous rock mass 
located below. Further, four horizontal strips are distinguished by the line 
z= 0.5, z = 0.85, z = 1.5, and z = 3.0 kilometers. Now it is easy to code 
the geological schematic since we have already written up the location and 
dimensions of the 14 contacts (Table 17). 


e TABLE 17. 
249 2 ee ee 
8 ae Paha No of the perturbing body 
E a a “ ae Paras i po 
(oP era | | 8 9°} 10 i4 
aoe. | ea oe 


—0,5 ;—0,514+-0,5 }+4+0,5| 0,4 | 0,4 | 0,4 | 0,4 
| 0,5 | 0,25] 0,5 | 0,25) 0,5 | 0,85) 1,5 
0,85] 0,5 | 0,85] 0,5 | 0,85) 1,5 | 3,0 

2,1 | 2,6 | 2,6 | 4,75) 4,75) 4,75] 4,75 


Commas indicate decimal points. 


In curve Veg we will fix the most descriptive points. As in the previous 


example, they are not selected uniformly. Wherever the anomaly is monotonic 
and has a uniform gradient, the points are selected more rarely. The points 
were selected more thickly at the most salient points of the curve. In all, 20 
points were fixed. A list of those is shown in Table 18. In order to calculate 


Pein? we will consider that the selection error should not exceed 5 eoetvoes. 


In this case F,._ = 1000 eostvoes- 
fin 
TABLE 18. 
= To re a =o mes = 
No of the | ° Pee a eee (Ree 5 ee Gnesi 
| points ; 5 3 i i ae fee ‘ aun eee eine 
¥ 0, | 075] 1.0 | 1.25] 15 | 475) 20 | 2.25] 25 | 275 
| Vee | 17 | 38 | 54 | 82 | 121 | 1077} 40 3 | —5 —24 
{ 
ae aes " | Iz |» | 18 |: | | 18 | 19 | 2 | 
ai ; 
x 3,0 325 | 35 | 40 | 475] 50 [5.25 6,0 | 7,0 8,0 
Vxy2 | —80 | —86 | —35] 33 | 103 | 101 | 75 | 36 | i8 i 


Commas indicate decimal points. 


The results of solving the problem are shown in Figure 10. Here the anomaly 
ee is cited which has been stipulated by the geological model hypothesis (first 


approximation). 
49 
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Figure 10. An Example of the Selection 
of a Geological Cross-Section (The Two- 
Dimensional Case) According to an Anomaly 
in the Horizontal Gradient of Gravi- 
tational Force. 


We will illustrate the solution of the problem by one more selection of a 
fairly complex geological object. Figure lla shows the density cross-section 


which has been constructed on the basis of available hypotheses about the 
geological structure of the region. 


We will fix the start of the coordinates on the line of observation. We 
will direct the axis Ox along the line of observations, and we will direct the 
axis Oz vertically downward. The area where the perturbing masses are concen- 
trated is divided into 11 horizontal strips by the straight-line z =h.. In 


each strip the limits of the perturbing masses are indicated by the straight- 
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lines xX = oe The density values of the rocks are written inside. The 


boundaries of the masses are enumerated (in the figure this numeration has been 
omitted). In the first and second strips, the boundaries are grouped together 
by nines (No. 1 - 9 and No. 10 - 18), in the third strip they are grouped 
together by 13 borders, etc. In all, 77 contact boundaries have been fixed. 

66 points have been fixed for selection. One notes at once the rather large 
difference between the observed anomaly and the results of solving the direct 
problem from the initial model (especially in the area of 50 - 55 kilometers). 
This is evidence of the fact that rather significant errors were allowed during 
the construction of the initial variation of the geological model. 


Figure l1lb shows the selected cross-section. The theoretical and observed 
anomalies coincided quite well. Their lack of agreement in small areas can be 
explained by heterogeneities in the upper portion of the cross-section which 
were not taken into account. This example is a good illustration of the change 
in the geological model of the first approximation and its transformation into 
another model. Completely satisfactory agreement was attained between the 
observed and theoretically calculated anomalies by this process. 
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Figure 11. An Example of the Selection of a Geological 
Cross-Section According to the Gravitational Force 
Anomaly: a) Initial Approximation, b) Solution. 
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3. Determining the Contours of a Perturbing Body With Partially Known /54 


Parameters 


Additional requirements for searching for the contours of perturbing bodies 
often arise when solving practical problems. 


Some parts of the cross-section are determined quite exactly on the basis 
of geological data or the data of other geophysical methods. These parts of 
the cross-section should not be changed during the process of selecting the 
contours. This means that among m steps by which the geological model is 
described, there are those in which the parameter g. does not change during the 


computation process. We will divide all of the parameters of the geological 
model (steps) into two groups. We will include in the first of these groups 
the projections belong to this group. The second group will comprise the 
remaining projections, all the parameters of which are constant. Now formulas 
(3.1) and (3.1la), by which the observed anomalies are approximated, should be 
presented in this form: 


' Ag (x) =k p> Fp hy Hy dp ER 2, F(a; hy, Hy dy, x), (3.6) 
a 5 ees sagen oe 
Vaal) = 86H hi Hn dye) +k BOC) hy Hy di, 2. ara) 


Functions f and » are determined by relationships (3.1) and (3.la). In equality 
(3.6) and (3.6a), the second part is only the function of the coordinates. 


We will now turn to (3.2). The function can be written as 


Ta a ie 
Fm 3 [dgayicttd—k 1 > / - 


j=m,t+l 


a3 2 (rene —A 3, i - $i] | 


or /55 


es p> (Agopie(*” = Agineo Ol» (3.7) 


where now Ag (x) is determined by equality (3.1) and is a function of only 


theo 
those projections in which the parameter qd is a variable quantity. The function 


Agibse (x) is determined by the relationship 


Bee De Sees as eee 
= d, (3.8) 
» Agobse(*) a Agobsel*) —k Dd fon hy, Hj, d,, x) 
and, is essentially a residual anomaly. The gravitational effect of the geo- /56 


logical projections, the parameters of which are known and cannot change, are 
excluded from the observed anomaly. 
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Small additions to the previously described program make it possible to 
fully automate the calculation. The parameters of the geological model should 
be fed in by using two massifs. One includes the vectors among the components 
of which there are variable parameters, and the other includes the vectors in 
which the components are constant. 


The second massif of initial data makes it possible to compute the second 
term of equality (3.8) and then function Agr se (x). In the future this function 


will be used as Ag (x). Everything mentioned relates fully to function {57 


obse — 
(3.2a) where the anomaly in the horizontal gradient is examined. We will 
illustrate this presentation with an example. Let there be in the geological 
model which was examined in section 2 (Figure 10) a well-known position for the 
line along which rocks with a density of 2.9 make contact with their enclosing 
medium. Thus, during the selection there is no basis for changing the position 
of the line. We will fix the parameters which describe this contact, and we 
will solve the problem by selecting the contours again. In this case, Table 17 
should be replaced by Table 19. Here the projections, the position of which 
should be changed, are clearly indicated. The parameters which describe the 
location of the contact line are written up in a Separate group. These quan- 
tities should not be changed during the selection process. 


7 TABLE 19. 

es uU 

Oo a} ee ee Se : ae. a ie ole 
5 of perturbing body 

ira be es Se ea Pt aad 

5's fo osb, aes 

ag pl ke “ash 98 


Commas indicate decimal points. 


In order to solve the problem, 20 points have been chosen on the anomalous 
curve, as was done in the previous problem. This means that the data of 
Table 18 are included in the initial data, in addition to the data of Table 19. 
The results of the solution are shown in Figure 12. Figure 13 shows the solution 
of this same problem for the gravitational force anomaly. The position of the 
contact on the right in the geological cross-section has been confirmed. In 
the model of the first approximation, a group of rocks, the contours of which 
are to be found, has been given in a more tentative manner than in the previous 
problems. 
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Figure 12. An Example of Selecting a 
Geological Cross-Section on the Basis of 
an Anomaly of the Horizontal Gradient of 
Gravitational Force. 
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Figure 13. An Example of the Selection 
of a Geological Cross-Section For a 
Gravitational Force Anomaly. 


CHAPTER IV 


DETERMINING THE CONTOURS OF THREE-DIMENSIONAL PERTURBING 
MASSES BY USING GRAVITATIONAL ANOMALIES 


1. Formulas for the Computation 


For an observed gravitational field, it is necessary to determine the 
Spatial distribution of a perturbing object which is finite in space. 


Everything which has been stated earlier about determining the contours of 
a two-dimensional body can be generalized to the three-dimensional case. We 
will look for the form of the perturbing bodies as a selection of straight 
steps which are a in space. Each of these steps is described by 6 para- 


meters (o, h, H, Li, t 2 d) The values of these parameters have been shown in 
Figure 1, 


We will consider that the excess density is a constant quantity within each 
of these steps. By changing the density value from step to step in each case, 


it is possible to describe the configuration of a sufficiently complex geological 
body of variable density. 


The gravitational effect of the masses included in such a step can be 
obtained by solution of the direct problem for a parallelepiped. It is sufficient 
to presume that parameter xX, 7% 


If there are m steps, then the observed anomaly can be approximated in the 
following way: 


For the gravitational force anomaly 


ie = os pa = ee 
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For the anomaly of the horizontal gradient: 


oe (4.2) 


where 


In the observed field we will fix n points with coordinates xX; and Yi- 


We will structure the same function as (1.2). In order to determine the 
parameter a we will write the values of the derivatives: 


For the anomaly of the force of gravity 


an ah n WT T: 
Ky = = Dh 2, See, 4 py oat yd X 
: _ bol, ae 
as (yp = yd + APUUCys = 9) #454 
iy Wu ied aT — 7 pa 
a ee 


ta$ doay, 
ere 


ee f (4; = xpi es, (dj — x;)? es 


Ty [ar —uo + AAI a (Gy — 9 + APTA? 
bye de ae 
a aria aa My — yi) + AF AF 


45 | 
Hy ola Mee —0| ar —ar|+ (4.3) 


a ye (ai — yo ar +A 
a By = HP (las — yi? + H?(A?)2} AP 
LE (hy — yi) My — yd? + 1) 

y— x xl? (yi — yO? + Hi. TAP?) Al? 


FE (lay — vd lg — 4o* + 3) £60_ 


a eee 
a Ay = xi)) (ly =O: + i BAM ay A + 


ly — yd Mh = + Hil 
cr = x0 hy yg ny(Ay) AP 


57 


For the anomaly of the horizontal gradient 
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2. Selection Methods 


The methods of selecting three-dimensional bodies do not differ basically 
from those described for two-dimensional bodies. The first variation of the 
problem provides for a search for the contours of the perturbing body in the 
coordinate plane x0z while taking into account its dimensions in space. If 
some of the objects are not intersected by the coordinate plane, then a pro- 
jection of the contours is sought. 


We will use examples to illustrate all our methodological means. 
The Gravitational Force Anomaly 


Let a gravitational force anomaly (Figure 14) be given. As is evident, 
the perturbing object is nearly isometric, making it impossible to use the two- 
dimensional problem. We will fix the beginning of the coordinates in the center 
of the anomaly. We select the axis Ox across the course being noted. Then it 
is natural that the axis Oy should be directed along the course, and the axis 
Oz should be directed vertically downward. Figures 15 and 16 show the curve of 
the gravitational force anomaly. We will draw columns like the Barton measuring 
grids in planes x0z and yOz. Let it be established that the initial depth can 
be accepted as equal to ho = 0.25 kilometers, and the final depth is H = 4.0 /61 


kilometers. We will divide this whole region into 6 strips by using the planes 
Ae = 0.25; 0.5; 0.875; 1, 375; 2.0; 2.75; 4.0 kilometers. We will solve the 


problem sequentially. At first we will make approximation in the plane x0z, 
then in y0Oz, etc. 


For the zero approximation we will select the bar (parallelepiped) which /62 


is bounded by the planes x= -1.75 kilometers, xX, * 2.0 kilometers, y," 3.0 


kilometers, and ce +3.0 kilometers. 


In the given case, in planes x0z and yOz the number of unknown parameters 
(d) will be the same (12 unknowns in each case). 


Let us examine curve Ag (x). We will fix the most descriptive points on 
this curve. The coordinates and values of the anomaly are given in a table 
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(Table 20). Table 21 shows the geological model. Of every six parameters, the 
first five are constant, and the sixth dis variable. Tables 20 and 21 show 
the results of the solution. The calculated field agrees fully satisfactorily 
with the observed. 


Now we will move to selecting the 
contour along axis Oy (Figure 16). The 
parameters which describe the space 
have been taken from the previous solu- 
tion (along axis 0x). Tables 22 and 
23 show the coordinates of the points, 
the values of the anomaly, and the 
parameters of the geological model. 

The results of the solution are also 
shown. The first cycle of calculations 
is concluded with these calculations. 


Analyzing the results of the cal- 
culations, one can state that the com- 
putation was done with considerable 
Figure 14. The Field of Gravi- exactness. The selected and observed 
tational Force. anomalies agree completely satisfactorily 

with one another. However, the contours 
in plane x0z have been selected with 
extremely approximate values for the dimensions of the masses which are anomalous 
in space, As one can see from Table 21, the uppermost block was computed under 
the condition that its space ranges from -3.0 to +3.0 kilometers. The values 
of these parameters have been computed according to the profile along axis Oy 
and are -1.83 and +1.13 kilometers. In order to clarify how this is reflected /64 
in the results of the solution of the problem, it is necessary to repeat the 
computation, making a unique second cycle. Thus, in Table 21 parameters Ly and 


L, have been taken from the results of the computations of the first cycle along 
axis Oy (Table 23). 


The values of the parameters computed in the second cycle Ox are shown in 
Table 21. 


Now the parameters in plane y0z can also be defined more exactly. After 
changing parameters Ly and Los we will carry out a second cycle of computation 


along axis Oy. The results of these computations are shown in Table 23. When 
analyzing these data, it has been established that now the discrepancies in the 
desired parameters are small (they are less than 200 meters). 


A third cycle was also carried out. Its results differ little from the 
second and are not shown in tables. 


The results of the calculations are shown in Figures 15 and 16. 
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Figure 15. Selection of the Profile 
Along Axis Ox. 


Figure 16. Selection of the Profile 
Along Axis Oy. 


TABLE 22. 
eee ee EE a ee aS 
Coordinates - Value of va - 
a. 2-q in the initial Selected in 
_| Observed approxima: ' the first 
as +ftion (1 cycle)k cycle : 


. m 


2nd cycle 


*. 0,78 
“Eo | 0 476 
—45} 0% 3,53 
—30| 0 — 
og aa 10.45 
—15 0 Z 12,88 
—10 |. 0 14,33 
04° 15,04 | 
nel 8: 14,87), 
lol. 0° 14, bee 
20! oO 10,47 
30] 0 ee 
101 0 44st 
55 0 215. 
go} 0 0,75 


Commas indicate decimal points. 


TABLE 23. 


Parameters of 
the geological 
1 : 


) scheme boa 
ge a —0,25 
ys h 0,25 
.H 0.5 0,875 0,875 1,375 1,375 
oh —0.93 | —1,05 | —l,05 | —1,32 —1,32 
sie ls 0,65 | O46.t.° 0:44 | 0.68 0,68 
dassip ¢. 3.0. [sr3,0.° [ ~°3,0° | —3,0- 3,0 
4 


ee 
lcycle “}? -=}; 
II ‘cycle - f rly 


Parameters of 


“025 | —o25 | 025 | —0,25 
0.5 05 | - 0.875 | 0.875 


Selected in the 


the geological 8 9 | “to | it 12 
scheme et 
Pt 0.25 | 025} —025 | 0,25.) —o,25 
h 1,375 2,0 2,0 2,75 2,75 
_H 2,0 2,75 2,75 40 | 4,0 
lh 14 | —07 | —07 | —1ez | —187 
E 170-}' 228 | 228] 265 2°65 
d, ssig .. 3,0° —3,0 3,0 —3,0 3.0. 


ie ake 
dealcul 
I cycle 


II cycle 3,06 


—3,52 


Commas indicate decimal points. 
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Figure 17 shows anomalies of the horizontal gradient of the force of gravity 
for the same geological object. The selection method remains the same as before. 


In Tables 24 - 27 data are cited for the computation of the inverse problems 
for anomalies Vee and Mos The parallelpiped, the boundaries of which converge 


with the extreme points has been accepted as the initial approximation. At 

the very beginning, using the approximated data about the dimensions of the 
‘body, along axis Oy, we determine its contour in the section x0z (Table 24 and 
25). Then we move to the anomaly along axis Oy. The overall dimensions of the 
body have already been defined more exactly by the previous computation, and 


the results are shown in Tables 25 and 27. The first cycle of computation is 
concluded at this point. 


Number of 
the Point 
I cy 
Number of 


Commas indicate decimal points. 


TABLE 25. 
Parameters of | co: 
the Geological 1 | 2 | 3 | 4 f: coe | 6 | 7 | g | 9 | Jo uy 
Mode - m 4 ! 
= 
eC 1 |. 1 1 fal 1 i 1 7. 
er: o5~ Tt 05 0875 [0,875] 1.375 | 1,375 20 } 2.0 | 2.75 

bee 0,875 | 0,875 | 1.375 f.1:375| 20 |. 20 275 | 2.75 |, 40. 

h —U75 | “4,78 |—175 |—L7e |—1.75 fais | 1B | —175 | —hIS | — 
“tole a 1,75 1,75: 1,75 1,78 1,75 ° += 1,75 4,75 1,75 a 
5 dasels —1.0 10 |—Lo 10 |—1.0 1.0 —t0 1,0 | Beal, 

. g epee 
dealeul. ' = Ae ; 
“peycle | —1,0 0,54 | —1,10 0,98 | —1,5) 1,77} —1,69 L96 | —1,68 190 |, ae 
I cycle | —093| 0,47 | —1,08 0,96 | —1,37 1,36 | —1,63 213 | 73 | 213 mm 


Commas indicate decimal points. 
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The second cycle again begins by selecting the contours in plane x0z. 


By 


using the results of the computations along axis Oy, we define more precisely 


the body in space. 
shown in the same tables. 


Number of 
the Point 


Parameters of 
the Geological]: 
Model a 


oy 


I cycle —1,60 150 | —1,77 
Wicycle ‘| —1,67 1,66 | —1,74 


“v] Value of the anomaly Vv. 


2. 

3 

4 

5 

6 

7 

BL 

9 |—1,0 | + 67 | 109 

10 |} —0.75} 40 75 

11 | —0,5 25 47 
Commas 


TABLE 26. 


and 
Ou 
ie 
$5 Bel y 3 
ooo =I oF ra) 
~~ a e oO 
D & OVA 3 
ae] . a 
48 12 - 0 0 y —2* 
86 13 | 05 —25) —47! —397 
118 4 1,0 —67|—109] —86 
148 15 1,25 | —104} —153} —117 
177 16 15 | —173) —213} —151 
207 17 1,75 | —200 | —260 | —178 
4 172 18 | 2,0 |} —179| —219] —179 
; 109 19 2.5 | —140| —127} —140 
83 69 | 20 3,0 |—115} —80] —116 
57 44 |) 21 5,0 —52! —19} —25 
35 26 
indicate decimal points. 
TABLE 27. 
5 6 ? 8 9 
1 1 i 1 i 
0,875 | 0,875] 1,375 | 1,375 | 2,0 
1,375 | 1,375] 2,0 2,0 2.75 
—15 —15 j—1,7 —1,7 —1,68 
1,77 1,77 1,96 1,96 1,90 
—1,75 1,75 | —1,75 1,75 —1,78 
—2,00 2,00} —2,17 2,09 | —2,08 
—2,82 2,96{ —2,71 2.69 | —3,42 


Commas indicate decimal points. 


The computed anomalous function coincides sufficiently well with the 
observed function, and there is no need to repeat the iteration process. 


The results of the selection are shown in Figure 17. 


cf in the 2nd 


LF 


cycle — 


The results of the computations in the second cycle are 


Here two graphs are 


combined, and therefore the valués..of the anomalies in the first approximation 


are now shown in the tables. 


model for which all the computations were carried out. 
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The figure shows the contours of the theoretical 
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CHAPTER V /69 


INTERPRETING ANOMALIES COMPLICATED BY LOCAL BACKGROUND 
CONDITIONS 


l. Posing the Problem 


It is very important to know the level of a normal field when interpreting 
gravitational anomalies. For a variety of reasons, we very often have no 
absolute general knowledge about a gravitational force anomaly. These reasons 
can be inexact knowledge of the normal distribution of gravitational force and 
extremely approximated knowledge of the density parameters of the rocks which 
make up the intermediate layer. As a result of this the Buge correction factors 
and many other things can be very inexact. 


When solving a direct problem for certain geological conditions, various 
investigators can obtain different readings for the field level. This can be 
illustrated by a very simple example. 


Let there be in the cross-section a limitless plane-parallel layer in which 
contact is clearly distinguished between two types of rock which differ in 
density. Examining this cross-section from the side of the less dense rock, we 
obtain a general positive anomaly for the force of gravity. If the investigation 
of the region is carried out from the side of the dense rocks and if the normal 
field is selected at that point, then the anomalous effect of the contact will 
be located in the range of negative values for the force of gravity. The first 
and second anomalies differ from one another only by their constant components. 
Many similar examples can be cited. 


We will now turn to questions of interpretation. The machine method of 
selection examined by us consists of comparing the observed field and the 
results of solving the direct problem. It is evidently clear that in a number 
of cases these fields can differ not only in details, but can have an absolutely 
different level reading. No attention is paid to this during manual selection. 
A simple parallel transfer (often completely automatic) solves the entire question. 
The situation is different with the machine selection method. Imaginary per- 
turbing bodies will appear when the calculated body is combined with the observed 
body. 


We will turn again to an example. Figures 18 and 19 show the results of 
selecting the contours of geological bodies for a gravitational force anomaly. 
The entire process was carried out using the computer method described in 
Chapter III. In the lower parts of the figures the geological model is shown 
which was constructed by the investigator as the initial variation (the first 
kilometer of the cross-section has been excluded from the investigation). The 
direct problem for this geological structure on the level of the field differs 
by 15 mgl from the observed field. 
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Figure 18. An Example of Automatic Selection Without 
Taking Into Account a Difference in the Levels of the 
Observed and Calculated Field. 


Figure 18 (the upper cross-section) shows the results of the selection. 
In order to compensate for the 15 mgl, horizontal intermediate layers with excess 
densities of +0.27 and +0.17 have been inserted into the geological cross- 
Section. These masses have been formed by the movement of heavy masses to the 
left. The horizontal dimensions of the masses which had a density of +0.1 and 
+0.2 have been sharply increased. Although the observed and selected anomalies 
coincided quite well, it is obviously difficult to evaluate positively the 
results of the selection. 


A completely different picture was obtained after introducing into the 
observed anomaly correction factors for the field level (Figure 19). The 
observed anomaly and the results of solving the direct problem, it is true, 
differ considerably, but they coincide in their asymptotic parts. Now the result 
of the selection differs from the initial approximation only in a few details 
(although these details are very significant). 
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Figure 19. An Example of Automatic Selection Taking 
Into Account a Difference in the Levels of the Observed 
and Calculated Field. 


Thus arises the task in the observed field of distinguishing two components, 
one of which can describe the geological heterogeneities, and the other the 
general local background. In its mathematical formulation, the task will consist 
of the following. Let an anomalous field ee — | y) be assigned. (This can 


be any gravitational or magnetic anomaly). a will fix in this field n points 
with coordinates xX; and Yi- On the basis of all information about the geo- 


logical structure, taking into account the field being interpreted, we construct 
a geological model in such a way that the anomalous effect from it can be com- 


puted. Let this effect be Med (x, y). Using the function f (x, y), 


approximate the regional component of the fields as ve reg = £f (x, y). Now we 


will compare the two groups of the function V onee (x, y) and V ties (s, y) + 


f (x, y). We will construct this expression: 


F= -3 Vase (Xi yi) —Viheo (xi, yi) —f(Xp y,)/. (s . 1) 
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For a certain geological cross-section and at the assigned coordinates of the 
points, the variable parameters will be contained here only in the function 


viag = f (x, y). Expression (5.1) can be rewritten thus: 


= — 


P= SIA 9 — Fey yl (9:2) 


where ths 
Ale 4) = Vobse (Xe Yi) Viet y,). 

If one takes into account the fact that the observed anomaly is given by 

the table and that the theoretical anomaly can be computed at all fixed points, /73 


then it becomes clear that it is easy to compute the function A (x, y) at all 
fixed points and to write these values in tabular form. 


Now we will turn to the function f (x, y). Having given a certain form 
to this function, it is possible to find from the conditions of the minimum 
of (5.2) the parameters which describe it. Thus in the very first approximation 
it is most convenient to approximate the local component by using a linear 
function with the form 


"Ba, y= A+ Bx+Cy. 


One can almost always find a geological explanation for this function. In this 
case (5.2) can be written in the form 


F = D A(x, y,) — A— Bx, — Cy,). (5.3) 
ive] 


Now the task is to minimize function F (A, B, C). 
2. Minimization Methods 


We will examine two computer models for determining the parameters which 
describe a background field. 


Solving the Problem Using the Fast Descent Gradient Method 


It is possible to minimize function (5.3) and to determine the parameters 
of the local background using the fast descent method. The computation method 
does not differ from those used in the problems examined earlier in this work. 
The desired values are determined by the sequential approximations 


wo —-—S 


At A” — Fa: (5.4) 
Beth = BY =, | Fs 
ct) co a 


The function F and derivatives of it depend on the geometric forms by which the 
geological bodies are approximated. 
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Solving the Problem Using the Linear Algebra Method 


The problem described earlier can be formulated in the following form. 


The anomaly computed from the geological model hypothesis aes (x, y) and the 


regional component of the field Vyeg (x, y) should fully correspond to the 
observed field, i.e. 


Vopse(X» ¥) = Viheo(% #) + Vreg(* 9): 


If the regional component is approximated by the linear function A + Bx + Cy, 
then the last relationship can be written thus: 


Vobse(*, 9) = Mineo (sy) + A+ Be+Cy. 
For the fixed geological cross-section 


A+ Bx + Cy = Vase (4) Y) — Ving he ) = ALE I. 


We will fix n points with coordinates Xx. and Y; and we will construct 
the system 


“APxXB+y,C =A (Xp Ys (5.5) 
t= 1,2, eneirel ge EGS 


If n > 3, then a re-defined linear system is obtained which is written in 
matrix form as 


Ka=A, (5.6) 
where 
Ixy A A(x, 41) A, 
K = Ix, Ys, ; ee B ; A= : (Xo, Yo) ‘ 
IX,n - Attn id. LA, 


This system can be solved using the method of the smallest squares. The method 
itself relieves us of the necessity of investigating the compatability of the 
given system. For each variation of the geological schematic and the system 

of assigned points, we will obtain the optimum solution to the problem. 


Without citing computations, we will simply indicate that vector a is 
found from solving the system 


ge he KS 
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Here K' is the transposed matrix K. As is known, this relationship always 
leads to a fully determined system of exactly the same number of equations as 
we have unknowns [45]. 


For our purposes the matrix for deriving P = K'K can easily be determined. 


It is always symmetrical: 


Lob 3320 ed n ae DTT 
pests Xq eee Xn Xa Me f= Dx; Yi Dxiyi |. 


t te ALI x | Lae Ses Dsl 


Vector K'A is determined thus: 


Kil age k = DA, 
KA=|x, x4... %,/[ 2 |=] DA 


4 Yo wee Yn A Dads 
n en 


Thus this linear system of equations is to be solved: 


nA+ DxB+ SyC= DA 
Dx A+ DxiB + Sey = VA, 
Dy A+ Sry.B+ Dy = VyA,. 


(5.7) 


If the observations are given in the one profile y = 0, then this system 
leads to 


nA + ¥x,B = ZA,, (5.8) 
DxA + MGB = 3 x,A,. 


3. Determining the Regional Background of Local Symmetrical Gravitational 
Force Anomalies 


For local symmetrical anomalies where geological bodies have a branching 
structure, the form of each geological object can be approximated by a sphere. 
If there are in all m geological bodies, then the anomaly from the perturbing 
masses can be written in the following way: 


m (5.9) 
; = Vv = Mily : 
ABtheo™ ¥) tea + OP 


j=l 


It is most convenient to express linear values in kilometers. We will decide 
to consider excess mass in units of 10? T. For instance if Oe 1019 T, 
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then it is necessary to insert into the formula M = 120 units. In this scale 
the coefficient K = 6.67 and the gravitational force anomaly are expressed in 
milligals. 


The calculation diagram can be described in the following way: 


1. Initial data. In order to solve the problem, it is necessary to assign 
a diagram of the geological structure. Its parameters are entered into a 
special table where the mass of each body and the coordinates of its center 
of gravity are fixed. 


In the anomalous field n points are selected which can completely describe 
all the peculiarities of the observed field. The coordinates of these points 
and the value of Ag are entered into another table. Thus the initial data for 
solving the problem are reduced into two tables. 


2. For the point, the coordinates of which are assigned, the function 
Ag neo (x, y) is calculated according to formula (5.9). The results of the 


calculations are printed out. 
3. The function is computed 
Ay = A (Xp Yi) = AB pce (Xr Ya) — AB eg Xe Yi): 
4. The coefficients of equation system (5.7) or (5.8) are computed. 


5. The system of 3(2) linear equations are solved and the regional back- /76 
ground is determined. 


6. The linear portion of the observed field is isolated, and the following 
function is sent to printout 


Ne ct y) = Ags (x. 4) — A— Bx —Cy. 


Now it is possible to analyze the solution of the direct problem. For this 
purpose, the two functions printed out are compared: Age ase (x, y) and 
AB theo & ¥)- 


Then there follows a selection of parameters for the geological objects 
using the method explained in the previous chapters. 


If the initial approximation and the solution of the problem differ signi- 
ficantly, then the determination of background parameters can be repeated. 


We will illustrate the calculation method with the following example. Let 
a field with an anomaly in the force of gravity (Figure 20) be assigned where 
four local objects are distinguishable against a background of the general local 
field. We will approximate these objects using spheres. The parameters of the 
spheres are shown in Table 28. 


72 


Figure 20. A Gravitational Force 
Complicated by Local Influences. 


TABLE 28. In order to determine the background, /77 
we will select 31 points distributed in the 
most descriptive areas of the observed fields. 


o 

bo wn 

be ef We will list these points in Table 29. 

a pag a ie (a 

oa? ? = s After calculation has been completed, 

a function is set up for the local background 
1 1,07 a 0 1 (3 - 0.5 x - 0.5y) mgl. At all fixed points 
a on : = ; a calculation is made of the rest of the 
4 ,135 2 0 0,5 field: 


Ag bee(% y) = Agobse(* y) c= [3 — 0,5% — 0,5y). 
Commas indicate decimal points. — - 
Only the geological heterogenities can now 
explain the observed anomaly. 


4, Determining the Regional Background of a Gravitational Force Anomaly by 


Block Structure 


We approximate the geological objects by selecting their straight projections 
which are finite in space. If m projections are given in the model, then the 
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gravitational force anomaly from the geological objects is approximated by 
formula (5.1). For two-dimensional purposes, formula (3.1) is used. The cal- 
culation diagram does not differ from the one described above. 


TABLE 29. 


‘ hobs og t : o Ju a 

oe e f= Pa | a 190w a ow a 

Sele TY 8 eee | Bd eS ae] wep ee 

g3|"| 7 | Pek}. | = [ge 3 
pos a a ee 

| 

| 1} 90 O | 10.73 12 j—1.0) —1 5,67 22) 2,5 )—1 5,13 

21.15 0 | 5,93} 13 }—0,6} —2 5.294 2 0 | 7,50 

3] 25 O | 4,57} 14 |—1,0 I 4,67 24) —1 0 | 6,32 

4 lay O.} 1,23] 15 O} 25) 2,78) 25] 2,5 |—2.5 | 7,13 

5|;u0 |—! 6,69) 16} 0.5 15] 4,35) 26) 2 |—2.5 | 8,66 

6} 1 —! 6,08] 17] 1,5 1 5.15] 27] 3 2 | 3,30 

' 7135 }—-1 2,94) 18 2 2 8.57] 28] 2 3 | 3,23 

' 8h) —2 6,64]| 19 3} 0.5/ 283] 29] 0 1 5,69 

g'2 —2 |10,57|| 20 4] 4 |—0,67] 30] 1 O | 6,32 

10°52 | —3,5; 5.42 2 | 6.06] 31] 4.5 |—4 3,00 
HW) --3 0 4,82 | 


=. 


Commas indicate decimal points. 


The calculation method can be illustrated by the following example. Figure 
21 presents the geological model of the region and the gravitational force 
anomaly. Behind peak No. 20 in the rock mass, which has a density of 2.66, 
there is a strip of very heavy rock. The task is to divide the observed anomaly 
in the area from the zero point to the peak at point 10 into two components: /78 


- Agobse(*) =; Agopie(*) +At Bx. 


In order to solve this problem, we will locate the descriptive points on 
a curve. Twenty of them are delineated (in our case n = 20). A list of these 
points is presented in Table 30. 


TABLE 30. 
38 
eS 2. 1 2 3 4 5 6 7 8 9 10 
x 0,5 1,0 2,0 2,5 3,0 3,5 4,0 4,75 | 5,25 | 5,75 
Ag 140 | 2,15 | 3,00 | 3,45 | 360 | 4,10 | 2,60 | 1,60 | 245 | 1,75 


Commas indicate decimal points. 
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The geological cross-section is represented by a selection of projections. 
In this case 59 projections are given in the cross-section. The parameters 
describing the location and dimensions of these projections are shown in Table 31. 
TABLE 31. /79 


Number of the body 


“& | +005 | —0,05 | +0.07 | —0,07 | +0,05 | —0,05] —0.04 | +0,04 
rm] 0,05 0,05 0,05 0,05 0,05 | 0,05] 0,05] 0,05 
{025 0,25 0,25 0,25 0,25 | 0,25} 0,25] 0,25 

i. ae 3,8 48 5,4 625 | 7,2 | 815] 8,25 
lag Number of the body ma 
| Bb & 

2. AS 9 10 a 12 13 14 15 
(1: ens Bees See Cane eee neeeeneS UPENeeN 
s.6° { 4007 | —007 | +007 | —007 | +022 | —0,19 | +0.95 

A 0.05 |- 0.05 0,05 0,05 0,05 0,05 0,25 

H 0.25 |. 0,25 0,25 0,25 0,25 0,25 05 | 

d 8,75 9,05 9,45 9,95 | 18,7 19,45 1,25> 


ee 


3 Number of the body 7 Fr 
ELE ae) ea : ie 
aQaey 16 | 17 | 18 19 20 21 22 3, | 
4H SO H ; a 
ge eS ee Ne 
0 0.05 | +0,07 | —9.97 | 10.05 | —0,05 | —0.04 0a! }-0,07 
A 0,25 0,25 0,25 0,25 0,25 0,25 0,25 0,25 
H 0,5 0.5 0,5 0,5 0,5 0,5 0,5 0,5 
d 3,6 4,6 5,2 6,25 7,15 8,15 8,25 9,0 
5} | Number of the body 
Soeur es ee, Wires er 
gana} 2 | Le | 26 | 27 28 | 29 | 0 
ds 5 0 : 
Log —0.07 | +007} —007 | +007 | —0,07 | +-0,07 | —0,07 
| oh 0.25 | 025} 025 | 025 | 025 | 0,254) 0,25 
i H 0,5 0,5 0,5 0.5 0,5 0,5 0,5 
1 d 9,1 9,45 9,9 11,6 12,15 14,0 16,5 
! 


Commas indicate decimal points. 
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TABLE 31 (continued) 


—; Number of the body 


v 


31 32 x Hu 35 36 37 


Parameter 


+022 | —0,19- 7-40.95" | 0.08'}..4+0,07 | 0,07] +£0.05 f-0.05" 


0,25 og6-oa-} 05 | 05 | 05 | 05 fps: 
0.5 O58} "LO2} Lo 1,0 10} 10] 40.. 
18,7 19,45 1,15 3,35 45 | 805) 63 | 71 
% : aN 
ae Number of the body - 
S| 38 40 4 a2 a 4 “1 45 
g . 


-—0,04 | +0,04 | +0,07 | —o,07 | +0,07 | —0,07 | +0,03 
05 0,5 0,5 05 0,5 0,5 0,5 
” 1,0 1,0 1,0 1,0 1,0 1,0 1,0 
8,15 8,25 9,05 9,15 9,35 985 | 19.0 


8 Number of the body 
o 
ESS 
Eo 46 47 48 49 50 Bl 52 53 
Bwopg 
a +0,05 —0,05 +0,07 —0,07 —0,04 | +0,04] +0,03 | +0,05 
A 1,0 1,0 1,0 1,0 1,0 1,0 1,0 1,75 
H 1,75 1,75 1,75 1,75 1,75 1,75 1,75| 3,25 
d 1,05 3,15 4,25 4,85 8,15 8,25] 18,75) 0,75 
6 Number of the body 
a 
Q 1 
5 2 ae 5 65 56 57 58 59 
£3 38 
a a i aad Pay a es ics oe Ge ra 
: 0] —0,05 +0,07 —0,07 —0,04 -+0,04 +0,03 
+ oh 1,75 1,75 1,75 1,75 1,75 1,75 
Hq” 3,25 3,25 3,25 3,25 3,25 3,25 
d 3,0 4,05 4,6 8,15 8,25 18,0 
Commas indicate decimal points. 
The results of the calculation are as follows: A = - 0.024 mgl, and 


B = 0.0665 mgl/km. Figure 21 shows the observed anomaly and the anomaly from 
which the linear components of the field were excluded. 
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Figure 21. An Example of Excluding 
the Linear Component of a Gravitational 
Field. 
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CHAPTER VI /81 


DETERMINING THE PARAMETERS OF A PLANAR CROSS-SECTION 


1. Posing the Problem 


Earlier we examined problems in which we searched for positions of geo- 
logical objects at which the difference between the observed and the theoretically 
calculated anomalies would be the least possible. Some parameters were considered 
constant. This relates primarily to planes. Let us look at Figure 22. Here 
an observed anomaly of gravitational force is shown. A diagram of the geological 
structure (cross-section along the profile) has been constructed on the basis 
of the totality of geological data, while taking into account the anomalous 
field. The direct problem has been solved for the given distribution of per- 
turbing masses. We will compare the observed anomaly and the solution of the 
direct problem. Even the most cursory analysis shows that a single shift in 
contours cannot give us a satisfactory result. It is certain that, analogous 
to the example examined in the Fifth Chapter, we can obtain here a solution 
which obviously is far removed from the basic geological diagram. In addition, 
it is necessary to change the rock densities only slightly for the observed [82 
and calculated anomalies to coincide completely satisfactorily. After changing 
the planar parameters, it is possible to move on to selecting the contours of 
geological bodies. Thus arises the task of minimizing function (1.2) where 
all the parameters in the geological diagram except for the density have been 
given and are considered known. 


agfirst approximation 
a 


Figure 22. Example of the Selection 
of Planar Parameters. 
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In its general form, the task can be formulated this way. Let a geo- 
logical diagram be constructed on the basis of data about the geological 
structure of the research regions, and with allowance made for the observed 
gravitational field. In this process the gravitational effects from it can be 
represented in analytic form: 


na 


AB eolt Y) = Bil (hs Ys Pate Pats = + Paid 


(6.1) 


All further statements will be true for any other gravitational or magnetic 
anomaly for which the additive property is valid. In this case the geological 


diagram has been approximated by m elementary bodies (for instance, projections). 


The parameter a (k = 1, 2, ..., s) characterizes the location and dimensions 


of these bodies. In order to solve the problem, we will use the method we 
examined earlier. Then this functional is subject to minimization: 


obse eo 


na ia 5 he o 

F — bs [Ag (xp 4.) — Ag. (Xj, CA (6.2) 
In order to set up this functional, n points have been fixed in the field 

of observed values 2 (x, y). These are descriptive points where the field 


has extreme values, turning points, gradient changes, etc. 


Let us turn to (6.2). If one takes into account the fact that points 
(x; 5 y,) are given, Ag neo (x, y) is expressed by formula (6.1), and all its 


parameters except as are given, then the functional F depends only on the 


densities. Representing them as components of some vector {o.}, one can write 


tad 


TE SENG Gesonen As (6.2a) 


The expression (6.2a) can be minimized according to the computation diagram 
examined in the first chapter. We will write only the expressions for the 
derivatives 


Fo = —2 SA speel%ir Yi) — 
9% 2 Peobse (6.3) 
— Ag theo ti yd) - F(x Yi Pits Paty +++» Psi) 


In connection with the fact that, by their nature, planar parameters for 
homogeneous elementary bodies are entered by means of a linear multiplier, 
another approach to minimization (6.2a) is possible. In this process it is 
necessary to somewhat change the way the problem is posed. Having fixed n 
points with coordinates xX; and y;» we can compare the values of the observed 


and calculated anomalies from the geological diagram 


| Ag obselXi- ¥;) =) Of (Xi Yes Pit Pofr +++» Psi) = Dail; (X;. Yj). (6.4) 
| @t ieee Il isl 
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Expression (6.4) is a system of n linear equations with m unknown parameters 
o.. It can turn out that nis greater than m, when we would be dealing with a 


re-defined system of linear equations. If nis less than m, then (6.4) will 

be a system without supplementary definition. Both in the first and in the 
second cases, the system can be solved by using the method of the least squares. 
We will write it in matrix form 


, Koma (6.5) 
where A(t 4) fem) - ; © tin (Ms ) 
hi (Xq; Y2) i, (*2, Yo) eee bin (Xq Yo) 


es © © © © © 8 @© 8 oe « 


O71 Agobse*t y:) a 
@ 5, . eS Ag bse (tes ¥;) = - : 
a Agobse(Xn Yn) Ag, 


We find vector o from the fully determined system of linear equations 
K'Ko = K’a, (6.6) 
where K' is the transposed matrix K. 


As is evident from what has been said, the calculation process in its detail 
depends on which elements are used to approximate the geological model. The 
perturbing bodies are broken down into a sum of straight projections. Here two 
cases are distinguished, i.e., two-dimensional and three-dimensional geological 
models. 


The Two-Dimensional Case 


The geological cross-section is given by selection of straight projections, 
i.e., four-dimensional vectors with coordinates to,, hy, Hie ds} () she 2s ehes 


m). For a gravitational force anomaly the form of the function f (x, hy, He ary, 


can be determined from formula (3.1). If the calculations are carried out for 
the Ve horizontal gradient anomaly, then the form of the same function will 


be established from the formula (3.1la). 
The Three-Dimensional Case 
The geological model is given by selection of straight projections which 


are finite in space, i.e., six-dimensional vectors {o., h., H., 2,., 2,., d.) 
Gea. 2 i DD. ite ea A el 
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> 


For a gravitational force anomaly the form of function f (x, y, No a 


aee OTE d5) can be determined from relationship (4.1). If the calculations 


are done for anomaly Vege then the same function can easily be written from 


formula (4.2). 


2. Examples of Solving the Problem 


We will illustrate the solution to the problem by the following example. 
Figure 22 shows a geological model constructed on the basis of a fairly detailed 
study of available geological material and an analysis of the observed field. 

In the region four geological objects with densities of 2.65; 3.05; 2.87; and 


2.76 g/em” have been distinguished. The density of the enclosing rocks is 


accepted as being equal to 2.71 g/cm’. The dimensions of each body in space 

are determined on the basis of a map of the anomalous field. The distribution 
of the force of gravity was studied on the basis of profiles which were selected 
over the extent of each object. The configurations of the geological bodies 
located in the research area are represented as a whole by a sum of steps which 
are finite in space. In all, 84 steps have been given. 


Forty-one points were used to interpret the gravitational force anomaly in 
the profile. The direct problem has been calculated for the given model. The 
selection of contours of the perturbing bodies was made according to the method 
described in the fourth chapter. The results of these calculations are shown 
in Figure 23. The outlines of the perturbing bodies changed, especially that 


of the rocks with a density of 2.87 g/cm, The horizontal thickness in the 
section near the surface does not agree with known data. 


49 first approximation 
F enel 
/ ee 


_ 
/ 4gobse ™ 


Figure 23. Example of Selecting Contours 
Without Changing Density Parameters. 
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We will compare the observed field and the result of solving the direct 
problem. They basically differ in amplitude, but the left and right asymtotes 
of the curves are relatively close. This makes it possible to propose that the 
basic errors in the geological model are caused by density properties. 


We will select density characteristics. The new values of the excess den- 
sities are shown in Figure 22. As is evident, some of the geological objects 
are not homogeneous. Thus, in one of the bodies a gradual increase in density 
is clearly visible the more the depth increases, and in the others individual 
heterogenous rocks are distinguished. By itself, this material required a 
geological interpretation. In the same figure (Figure 22) the selected anomaly 
is shown. It is certain that it differed in detail from the observed anomaly. 
It is fully understandable that we could not obtain a good agreement between 
the observed and theoretical curve; after all in this case only individual 
densities were selected. 


Now the results of the calculation will be subjected to analysis. The 
researcher can restructure the model somewhat, define its parametric value more 
precisely, and then move on to determining the configurations of the geological 
bodies according to the method examined earlier. Having assigned the values of 
the excess densities, we will again select the configuration of the geological 
objects. The results of the calculations are shown in Figure 24. 


49 selected } 


Figure 24. An Example of Selecting Contours 
With Allowance Made for New Density Parameter 
Values. 
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CHAPTER VII /86 


AN AUTOMATED SYSTEM FOR INTERPRETING GRAVITATIONAL 
ANOMALIES (THE MINIMIZATION METHOD) 


1. General Remarks 


As is well-known, the inverse problems in gravimetric survey, in their 
most general form, are associated with incorrect problems in mathematical 
physics. There are two approaches to solving these problems. One of them is 
based on special regularizing methods [23, 43, 62-67], while the other approach 
consists of pre-selecting the class of the solution. All the parameters of 
the geological features are found in the selected class. Then they can be 
determined in a uniform manner [62]. 


Use of the functional minimization method to calculate the elements of the 
occurrence of geological features ensures adherence to the specific class of 
the pre-selected problem. This is the only way to regularize the incorrect 
problem. 


The fast descent method used to solve inverse problems is similar, while 
the calculation method described in the first chapter works only for monotonic 
convergence. 


If, when interpreting a specific anomaly, we proceed from various hypotheses 
about the geological structure, i.e., in each case the initial geological model 
is different, then the results of the minimization will naturally also be 
different. 


It would evidently be erroneous to assume that different geological 
hypotheses for the same anomalous field would lead to different minima of one 
functional, and that one of these would be global. The fact is that in each 
case a different function is subject to minimization. In form, these are all 
the same as (1.2), but in content they depend on the constant parameters not 
subject to change. 


We would like to emphasize that often an unsuccessful solution obtained 
due to a poor initial approximation will channel the interpreter's thoughts 
into a new area. In this case a new model of the geological structure is 
constructed, and one begins again to solve the problem. 


2. One Possible Means of Selecting Initial Approximations 
When solving the problem, it is necessary to minimize function F = F (Py3 
Po3 ieee PY which depends on m parameters. In order to solve the problem, it 
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is necessary to select the vector ee (j = 1; 2; ...3; m). The value of this 


vector's parameters can be Saiculated in the following way. For each of the 


parameters, let its upper and lower limits of possible values (pi dy and {p, 


be determined. This can be done without any kind of res ena. 
(0) 
jH 
points in m-dimensional anace for the parameters. The straight line L which 
unites these two points should be drawn through the entire area and be diagonal 
to it. Therefore, it is natural to search for the points of the null approxi- 
Mation of this straight line. Each point on the straight line has its corre- 
sponding value F. It is necessary to select as the point of the null approxi- 


mation that where F = Fin’ In order to find these points, it is necessary to 


solve the equation ~> = 0 or, in its expanded form 
q Eva p 


The totalities {p.°} and (PSs dy can be regarded as the coordinates of two 


ap j oF oF = 
. a -_ oe ata a Ope COS @, Sates on Foy COS Am = 0. 
Here G)s Gos +++, A are the direction cosines of straight line L which are 
expressed by 
(0) (Uy —— 
cos &; = Fin Fin Pin (j = 1,2, ..., mm) 


Function F along straight line L can be expressed as the function of one para- 


meter 2, and for this purpose it is sufficient to establish Ps = ag + Z cos ae 


Now it is possible to find the function's minimum by using the ordinary method. 


The sample method may be used to look for the roots of function f (2) = a 7 
0. This method, which is well known in mathematical analysis, is based on 
the fact that the interval (0, L) which is included between the initial point 


Qo) and.the terminal point (pig: is divided into s sections (a,, b,) where 


a =0, b =L, a. =b. .. 
fo) Ss 1 i-l 


We will find the intervals where f (a)*f (b) < 0. We will calculate (2 5 b). 

+ a+b 

2 2 

is the root of the function. The calculations are finished for this section. 

a+b 
2 


> 0 we examine 


If the numerical value of |f 


| is sufficiently small, then Ly = 


The search for the following sections continues. Otherwise, when f (a) f 
< 0, we examine the interval fa, 2 : e _and when f (a) f 2 5 5 
a+b 


5 » b}. It is easy to apply this method on the computer. 


the interval 
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When all of the roots of the equation have been found, it is possible to 
calculate the initial values of the parameters 


= a 
WPidinitian 0 J T text COS Qj. 


It is obvious that there will be as many initial approximations as there are /88 
roots of function F' (2). For each root a value F (Z) should be found. The 
most probable values are those at which F (Z) has its smallest value. However, 


it should be regarded as worthwhile to search for a solution at all values bss 


It is possible to accept as a final variation the one which best satisfies the 
known data on the geological structure of the region. 


If the analytic expression of function F is sufficiently complex, then its 
minimum can be determined by normal tabulation. In order to check the behavior 
of function F or straight line L and to find the minima points, we will divide 
a section of the straight line into s equal intervals such that they are 
sufficiently small that we are confident that the changes in F will be monotonic 
within the limits of each section. Now it remains for us to calculate the value 
of function F at the sequential points on the straight line, the coordinates 
of which are determined thus: 


ipl} = (pj best (Ap;}, 


0) 0) : 
eee. D; —p 2 3 ° . 
where App=—“——" , k is the number of the point, {Pi} = (pin). The minima 


= 


points are determined from the condition 


Fain > Fa < Feu. 


Each of these points, and there can be several of them on the same line, is 
sampled as an initial approximation. 


3. The Basic Components of the Automated System 


On the basis of what has been said, one can state that in the practice of 
gravimetric survey work an automated system of processing and interpreting the 
observed data can be constructed. This system consists of two independent 
parts. The first part combines the calculations for processing the observation, 
introducing various correction factors and reduction factors so the anomalous 
effect can be regarded as caused only by geological heterogeneities. Calculations 
of field transformation are also associated. The charts of the anomalous field 
and the various transformations are the final result. 


The first part of the automated system has been developed and is being 
used by many scientific research and industrial organizations. 


The second part consists of several sequential stages. It is intended to 
begin operation after a geological hypothesis has been constructed in the form 
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of sections of structural models, based on study of the anomalous fields and 


all information about the geological structure of the region being investigated. 


In this process a special role is attributed to the planar characteristics of 
all the various types of rocks included in the model. The preliminary results 
of the quantative calculations can be used to create a diagram of the initial 
approximation. One condition is introduced: one must know how to calculate 
the gravitational effect from the component diagram. For this purpose, the 
geological structure can be approximated by selecting cylindrical, spherical, 
or other bodies. In a complex geological situation the perturbing bodies can 
be approximated by selection of contacts (gravitational steps). For three- 
dimensional bodies the elements of the geological system can be a projection 
finite in space. All the components of the automated systems are shown in the 
block diagram. The initial data are two groups of information. The first 
includes information about the observed fields, and the second contains infor- 
mation about geological structure (the hypothesis is the first approximation). 


me 


te 


; f ‘initial Data: The Anomalous Field, 
~. [Geological 6 Moder ay ssl eiclatislene 


Determining 
| Background 
Function 


aoe + . 

¥ er 

Selection Selection of Contours 

~| for the Geological 
of Planes |< 
Model 

$$ ea aes 
: obese ee v= 
; 4 he 
| Analysis | cf 


Diagram of 
Geological | a 
Structure _} 


Block Diagram of an Automated System for Inter- 
preting Gravitational Anomalies (The Minimization 
Method). 


The following functional is subject to minimization 


Bee ke, yi Vine Me Yi) — Pe Mi 1, (7.1) 


At the very beginning Vites (x, y) is calculated, i.e., the direct problem is 


solved. This intermediate result has an independent value. Comparing the 
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observed and the calculated fields, the interpreter can introduce into the in- 
itial model additions and changes in case there are very large errors. This 
must be done in such a way that the succeeding stages are directed at auto- 
matically solving the problems which can be formulated this way: change the 
geological model in such a way that the observed and theoretically calculated 
anomalies coincide as well as possible. It is completely natural that, if 
there are insufficient data in the geological model, then it is impossible to 
obtain an effective solution to this problem. We will clarify this with an 
example. Figure 25 shows a gravitational force anomaly and the geological 
cross-section of the initial approximation. The investigator considers that 
the negative gravitational er fect is caused by a large block of rock with an 
excess density of 0.03 g/cm), The geological cross-section is approximated by 
several projections. Comparing the observed and the theoretical curve, we see 
that a geological body of an increased density has been omitted in the model. 
It should be introduced, and then one should more precisely define the initial 
approximation. This same figure shows the result of the automatic selection 
(the fourth stage of the system) without further definition of the primary 
model. It is evident from the configurations of the object that there was an 
attempt to introduce heavy masses within the body. It is impossible to draw 
the contours with confidence in this case, since the necessary information is 
lacking. 


‘ In the second stage of the automated 

system, the local components of the field are 
determined (if this is required). In formula 
(7.1) all values except the parameters of the 
background function f (x5, y;) have been 


\ Pre, 


\ “Sselect! assigned. It is almost always possible to 
\29 prima z approximate the local background by using a 
“ta pproximate / linear function with the form A + Bx + Cy. 

———— If there are geological premises, then the 


form of function f (x, y) can be made more 

complex. No matter what function we use to 
approximate the local effect, it is always 

necessary to thoroughly analyze the calculation 
results and to give it a geological explanation. 

In the opposite case some function can be 

selected as the background function, but its /91 
nature will be completely different. Let us 

examine this using an example. 


10 _Zy Wye 
WZ OOS LE 


We will assume that nfo contact has been 
introduced between two blocks of rock, i.e., 
a gravitational step has been omitted. Its 
effect is included in the observed field, but 
Perturbing Body Without In- not in the theoretical function. This means 
troducing Data About a Heavy that some function which complements this 
Gravitational Body Located omission will be determined. As a result of 
Within the Massif. the analysis, the interpreter should not make 

this mistake. 


Figure 25. An Example of 
Selecting the Contours of a 
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Now the local component is excluded from the observed function. The in- 
trepreter is required to solve the problem of the next step in the calculation. 
If the differences between the observed and the theoretical anomalies are 
basically amplitudinal, and if the planar parameters in the geological model 
have approximated values, then naturally the third stage should be the selection 
of the planes. The fourth stage will be the selection of the geological model's 
contours. If the interpreter has decided that the basic differences between 
the observed field and the calculated field can be explained by an incorrect 
configuration in the geological model, then the third stage is the selection 
of contours. Then follows the selection of the planes (if this stage is re- 
quired). 


Finally, the interpreter should analyze all the computation stages. Some 
stages can be repeated, and their sequence can be changed. If the geological 
model has undergone significant changes, then it can become necessary to re- 
calculate the background function; after all it was determined for a particular 
geological structure. Then the calculations can be repeated, i.e., one can do 
more literations. After each stage the divergence between the iterations is 
determined, and the interpreter decides whether to continue the calculations. 


4. The Results of Sampling the Effectiveness of the Automated System for 
Interpreting Gravitational Observations When Solving Geological Problems* 


The minimization method has been used with success to calculate the block 
structure of deep sections and surface areas of the Earth's crust, thus de- 
fining more precisely the structural peculiarities of the outline of the 
crystalline bedrock concealed under a thick covering of sedimentary formations. 
This method also makes it possible to study more in detail the deep differentiated 
intrusive massifs, and the quantitative characteristics of the formations with 
their complex sedimentary structure, these formations being characteristic of 
the geosynclinical and flexible parts of the Earth's core, for map making in 
difficult regions, and for solving other geological problems. 


In the following sections we will cite several examples of determining the 
quantitative characteristics of geological structures. 


Studying the Structural Forms of the Outline of the Crystalline Bedrock /92 


The territory on the western incline of the Ukrainian shelf was studied. 
Structurally, the region is characterized by a two-layer structure. The lower 
layer is represented by a multifolded complex of Precambrian crystalline rock, 
among which granitodes of the granodiorite type predominate. The upper layer 
is characterized by a horizontal or very gently sloping occurrence of sedimentary 
formations of loess (average density 2.12 g/cm?) and volcanogenic sedimentary 
rocks of the Volhynian Series (2.8 g/cm3). The characteristic features of the 
gravitational force field which corresponds to the given territory is the pre- 
sence of unique zone or belt nearly latitudinal in direction and consisting of a 


*This section was written in collaboration with V. A. Rzhamitfyniy. 
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series of relatively isometric positive anomalies of varying intensity. This 
zone does not appear in a magnetic field. The results of interpreting the 
data of the resistivity surveying indicate variations in the surface of the 
crystalline foundation. 


Figure 26 shows the geological cross-section of one of the areas in the 
anomalous zone. The crystalline rocks were investigated to a depth of 30 to 
65 meters by drill-holes bored into the center of the anomaly. In the marginal 
zone the resistivity surveying establishes this type of rock at a depth of 300 - 
350 meters. The presence of a strong gravitational force anomaly and a signi- 
ficant fixed value for the density gradient at the boundary between the crystal- 
line and the sedimentary rock makes it possible to determine the parameters of 
the perturbing feature. In the given case these are the outlines of a pro- 
jection of the crystalline foundation. The density characteristic of the rock 
making up the projection were determined by approximation, as was also the 
average density value of the rocks opened by the drill-holes (2.82 g/cm3). The 
approximation model was constructed by using all geological and geophysical data 
and was presented in the form of a geological cross-section in which the rocks' 
main characteristic was their density. The body which caused the anomaly’ was 
considered to be a three-dimensional object. The isometric form of the anomaly 
in the map served as the basis for this premise. The direct problem is solved 
from the model of the first approximation. The anomalous effect obtained in 
this process was compared with the observed gravitational force anomaly. There 
exists a considerable disagreement between the anomalies both at their maximum 
points as well as in the sloping regions. However, the disagreement is very 
regular, i.e., at every point the calculated anomaly remains below the observed 
anomaly, i.e., is small in absolute value. This type of discrepancy may be 
caused by the following factors: either the density of the rocks which make 
up the projection has decreased or the depth of its occurrence has not been 
correctly determined. Since the distance to the surface of the foundation was 


determined on the basis of drilling data and the results of resistivity surveying, 


it is possible to change the parameters of the projection only within the given 
depth limits, from 30 to 350 meters. Since the greatest disparity between the 
calculated and observed anomalies occurred in the zone of the maximum, and since 
the average value of the rocks' excess densities in the foundation has been 
determined in approximation, it became necessary to begin solving the inverse 
problem with automatic selection of the excess densities. As a result of the 
solution, the excess density of the crystalline rock was found to be 0.74 g/cm 
as opposed to the 0.70 g/cm? adopted in the first approximation model. In this 
process an insignificant differentiation of the density parameters within the 
limits of + 0.01 g/cm? was established on the basis of depth. This can be ex- 
plained by the lithological peculiarities of the rocks in this region. 


3 


From Figure 26 it is evident that the selection of density parameters 
yielded a satisfactory agreement between the observed and calculated anomalies 
in the zone of the maximum. The next stage of the automated system makes it 
possible to solve the problem concerning the configuration of the geological 
feature. This amounted to automated distribution of the elementary projections 
within the limits of the assigned depths at the selected density values. As a 
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result of the selection, almost complete agreement was attained between the 
gravitational force anomaly calculated from the cross-section selected and 
the observed anomaly. 


Nef oN? PII 


Figure 26. An Example of Selecting Density Para- 
meters and Contours for the Topography of the 
Crystalline Foundation: 1) Observed Gravitational 
Force Anomaly; 2) Anomaly Calculated From the First 
Approximation Model; 3) Anomaly Calculated From 
Model With Selected Densities; 4) Calculated Anomaly 
After Selection of Contacts; 5) Sedimentary Deposits; 
6) Crystalline Rocks of Projection; 7) Relief of 
Projection in First Approximation Model; 8) Relief 
Selected by Computer; 9) Drill Holes. 


Among the geological results obtained by use of the minimization method, 
mention should be made of the differentiated density boundary in the sloping 
portion of the projection and the gentler slopes themselves in contrast to the 
sharp, abrupt scarps shown in the first approximation model. The results ob- 
tained can be used in future investigations into the nature of the foundation 
relief. 


Figure 27 shows an example of a more complex case of determination of 
parameters which give a quantitative estimate of the projection of a crystalline 
foundation. Dense basalt mantles of the Volhynian series which occur in the 
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form of a stratum with a thickness of approximately 70 meters at a depth of 


around 70 meters also take part in forming the observed gravitational anomaly. 


The development boundary between these two formations is clearly fixed in a 


Magnetic field. 


As in the first case, all information about the structure of the region 
has been used to construct the first approximation model. 
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Figure 27. Example of Exclusion of Linear 
Components of the Regional Background and 
Selection of Relief Outlines of a Crystalline 
Foundation: 1) Observed Gravitational Force 
Anomaly; 2) Observed Anomaly With Allowance 
Made for the Local Background; 3) Amomaly 
Calculated From the First Approximation Model; 
4) Anomaly Calculated From the Second 
Approximation Model With Allowance Made for 
the Local Background; 5) Anomaly Calculated 
From the Selected Diagram; 6) Crystalline Rock 
of the Formation; 7) Basalts; 8) Sedimentary 
Deposits; 9) Outline of Projection in the First 
Approximation Model; 10) Outline Selected by 
Computer. 


In solution of the 


direct problem a considerable divergence was noted between the calculated and 


the observed anomaly, a divergence reaching 1 mgl in the left part, 4 mgl in 
the central part, and 3 mgl in the right part of the cross-section. 


The 
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difference of 2 mgl may be ascribed to incorrect selection of the level. A 

2 mgl difference between the extreme right and left parts of the anomalies 
being compared (where the effect is felt least of all) can be explained as the 
influence of a large anomalous mass located far beyond the limits of the given 
cross-section. This most likely is the effect of an abrupt general movement 

in the surface of the foundation from right to left. This effect can be in- 
cluded in the general local background, which in the given area can be approxi- 
mated by the linear function A + Bx. In this process the anomalies being com- 
pared are automatically combined at one level. The divergence between them 

in the central or maximum part can perhaps be explained by insufficient 
elevation of the central part above the general level of the foundation, or 
perhaps by the density differential of the rocks which make up the projection 
in the sketch. Since we had no confirmation of the second premise, it remained 
possible to reduced the occurrence depth of the highest part of the uplift. 
This correction was reflected in the second approximation model. The bend in 
the calculated anomaly was caused by a lack of mass of increased density in 
the aperture between the basalt mantle and the left slope of the projection. 
This discrepancy in the anomalous field can be eliminated by increasing the 
excess mass by moving the projection within the assigned depth limits. 


The direct problem is again solved from the augmented geological model. 
This is done in order to estimate the effects of the corrections been made in 
the first model. General agreement is noted between the calculated and the 
observed anomaly. After the parameters of the local background had been deter- 
mined, the anomalies agreed quite well in configuration (and, what is most 
important, they agreed in the maximum part); this is evidence of the fact that 
the occurrence of the anomalous features taken into account were determined 
correctly. The configurations of the perturbing features were then determined. 
The results of the selections are shown in Figure 27. 


Let us examine one more example. Figure 28 clearly demonstrates the 
capability of an automated system in interpreting a more complex gravitational 
force anomaly. Since there was no material available on the geological features 
which caused the observed anomalies, then, by anology with areas which had 
been subjected to more extensive study, and which were similar to the territory 
under investigation, the following concept of the geological structure was 
evolved. The anomaly is the aggregate effect of a thick tectonic zone made 
up of cataclastic rocks and causing a broad minimum, as well as of a projection 
of foundation rock causing a less broad maximum. The anomaly was primarily /96 
complicated by the very intense background caused by the general sharp sub- 
sidence of the foundation. Figure 28 shows the anomaly already free of the 
linear component, which was allowed for automatically as in the previous example. 


A highly simplified model of the geological medium which served as the 
first approximation model for quantitative calculations was constructed with 
allowance made for the following circumstances. The depth of the foundation 
surface, according to resistivity surveying data, is 600 meters. The cataclastic 
rocks within the limits of the tectonic zones are on the average 0.1 g/cm? 
lighter than the massive rocks, and the massive rocks themselves of the foundation 
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0.7 g/cm? denser than the sedimentary formations. It is evident from Figure 28 
that the anomaly calculated from the first approximation model differs con- 
siderably in absolute value from the observed anomaly, although it copies it 

in its general outlines. The solution of the inverse problem by automatic 
selection of contact points not only has made it possible to achieve good 
agreement between the anomaly calculated from the selected model and the 
observed anomaly, but has also made it possible to draw certain geological /97 
conclusions. Although the premise which served as the basis for creating the 
model of the geological structure has not been confirmed by actual geological 
Material, it nevertheless gives us reason to state that the given complex 
gravitational force field anomaly consisting of a broad minimum and a relatively 
Narrow maximum may most likely be caused by a thick tectonic zone which tends 
toward large values for the field of gravity. The maximum in this case may 

be caused by block uplift over the destruction zone or along one of the feather 
joints. The expansion of the tectonic zone in the upper part of the cross- 
section may be ascribed to the presence of a thick linear crust deriving from 
weathering of cataclastic rocks. 
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Figure 28. Example of Selection of a Geological 
Cross-Section in Interpretation of a Gravitational 
Force Anomaly: 1) Observed Anomaly; 2) Anomaly 
Calculated From First Approximation Geological 
Model; 3) Anomaly Calculated From Selected Model; 
4) Crystalline Rocks; 5) Zone of the Catclasis; 

6) Sedimentary Deposits; 7) Outline of Geological 
First Approximation Model; 8) Outline of Geological 
Model Selected by Computer. 
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Possibilities of Studying Folded Structures 


A region was studied within the limits of a territory which was characte- 
rized structurally by a two-layer structure. Metamorphic rocks of the gneissic 
series, as well as granitoids, have taken part in the formation of the lower 
multifolded layer. The upper layer is represented by cenozoic sedimentary 
deposits lying horizontally on the waterworn precambrian surface. The basic 
structural units which characterize the precambrian formation in the territory 
being investigated are the dome-shaped branchyanticlines. Their cores are made 
up of biotitic porphyroblastic granites (with a density of 2.63-2.67 g/cm®). 

In the other gravitational force field anomalies they are manifested in the 
form of relative broad gently sloping minima. The dome-shaped structures 
mentioned are flanked by amphibolous biotitic medium-grained migmatites 
(2.71-2.72 g/cm), The development areas of the latter are characterized by 
quite intensive (up to 3.0 mgl) maxima in the gravitational force field. Some- 
times tectonic contacts are noted between the above-mentioned various types 

of rock. Of the tectonic dislocations of a local nature, the sublatitudinal 
ones are considered to be the most recent. They control the movements of the 
dike complex, i.e., diabases (2.8 g/cm?). Biotitic gneisses (2.66 g/cm?) play 
a significant role in the geological structure of the territory. Intense 
erosion of the cataclastic rocks has caused the rugged topography of the sur- 
face of the crystalline base. In this process significant depressions are 
traced which are filled by light Cenozoic carboniferous deposits as well as 
areas where the weather-worn crust has developed to a great depth. These 
depressions are clearly reflected in the gravitational force field in the form 
of broad minima. 


The presence of a differentiated gravitational force field within the 
limits of the territory the geological picture of which is characterized by the 
development of different types of rock differing quite sharply in density 
characteristics (the average density was determined on the basis of samples 
from natural outcrops and from the core samples of drill-holes) have allowed us 
to attempt to analyze the abyssal structure by using quantitative interpretation. 
A first approximation geological cross-section was constructed for this purpose. 
In order to create this cross-section, all available factual geological material 
and the results of qualitative interpretation of geophysical observations were 
used. The depth of occurrence of the proposed structures was chosen tentatively 
with allowance made for relative comparison of individual areas of the gravita- 
tional force field anomaly along the profile. The entire geological model was 
converted to a density cross-section. It served as the first link characterizing 
the correctness of the hypothesis regarding the structural interrelationships 
among the petrographic types of rock distinguished. With this purpose in mind, 
the gravitational force anomaly was calculated from the cross-section of the 
first approximation and was compared with the observed anomaly. The results 
of this comparison are shown in Figure 29 and indicate that the geological 
hypothesis may be accepted as a basis. However, certain additions to the first 
approximation model are still required. Without such correction, it is impossible 
even in the most tentative approximation to solve the problem of the structural 
interrelationships of the rocks. In particular, the insufficient amplitude of 
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the maximum portions of the calculated anomaly indicates that the depth of the 
main structures made up of dense amphibolous biotitic migmatites is obviously 
small and should perhaps be increased considerably, possibly even multiplied 
manifold. This premise has made it possible to increase the depth of the 
structures to 1.5 - 2 kilometers, i.e., almost fourfold in comparison to the 
initial cross-section. 


The presence of narrow local maxima within the limits of the broad minima 
which correspond to the depressions and wide cataclastic zones often accompanied 
by dikes of diabases has made it possible to assume the presence of thin dike- 
like bodies of dense rocks. Even with a small thickness (on the order of 100 - 
200 meters), these bodies could be reflected in the gravitational force field, 
Since they are of considerable size and should be detected in the form of relief 
projections within the limits of the depressions as younger and denser formations. 
In addition, as research has shown, the dike formations within the limits of 
the region are usually characterized by a branchlike arrangement of several 
thin bodies, and this makes it easier to detect them in a gravimetric survey. 


Thus, in the revised cross-section several thin dike bodies are distinguished 
which are similar to the one reflected in the cross-section of the first approxi- 
mation. 


The presence of two relatively intense local minima within the limits of 
the maximum which formed by the structure of the left half of the cross-section 
has made it possible to increase the development depth of the body of cataclastic 
biotitic gneisses reflected in the cross-section of the first approximation, 
and has also made it possible to assume the presence of a body of the same 
rock, but of smaller dimensions, in the zone of the minimum near the left edge 
of the cross-section. 


Within the limits of the tectonic contacts, which are characterized by /99 
zones of cataclasis, the form of their minima and their intensity have indicated 
a somewhat lesser thickness of the largest zone and at the same time a great 
depth of development of the cataclastic rock. The intense gradient in the ex- 
treme right part of the cross-section indicates a steeper occurrence of tectonic 
contact and a corresponding increase in the depth of the migmatite structure 
at this end of the cross-section. 


The depth of occurrence of the foundation is very small (10 - 60 meters), 
while the difference between the crystalline and the sedimentary rocks reaches 
0.92 g/cm3. Even insignificant variations in the topography of the foundation 
are sharply reflected in the gravitational force field. Therefore, determining 
refinement of the parameters of the depressions was a very important task. /100 
Solving this problem was made easier by the fact that in a number of cases the 
rock mass had been drilled through. The corresponding data were taken into 
account on introduction of additions and revisions into the initial model. 
Because of the small thickness of the sedimentary deposits, the topography of 
the crystalline foundation is shown on a scale which exceeds by a factor of 10 
the vertical scale of the crystalline rock structures. These and other revisions 
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are reflected in the second approximation model (Figure 29). In order to estimate 
the effect of the correction factors introduced, the direct problem was again 
solved on the basis of the geological second approximation model. Comparison 

of the observed and calculated anomalies revealed relatively good agreement 
between their configurations, especially in the extreme portions of the force 

of gravity curve. This indicates that the correction factors were proper. 


Figure 29, Geological First Approximation Model 
(a), Geological Model Modified After Solution of 
Direct Problem (b), and Geological Model Selected 
by Computer (c): 1) Observed Gravitational Force 
Anomaly; 2) Anomaly Calculated From Density Cross- 
Section; 3) Anomaly From the Calculated Cross- 
Section With Allowance Made for the Locak Background; 
4) Cenozoic Deposits; 5) Gabbro-Diabases; 6) Amphi- 
bolous Biotitic Migmitates; 7) Biotitic Granites; 
8) Biotitic Gneisses; 9) Cataclasites; 10) Tectonic 
Contacts. 


Along with the agreement mentioned above between the configurations of the 
anomalies being compared, a divergence is noted in their absolute value, one 
which increases gradually from left to right along the relief. In the extreme 
right portion this divergence reaches 1 mgl. It is impossible to explain this 
by the geological peculiarities of the rock, since the same granites with a 
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density of 2.63 g/cm? are detected both in the left and the right portions of 
the cross-section. Therefore, this divergence can only be explained by the 
influence of anomalous features which are either located far beyond the limits 
of the area or which occur at a greater depth. This effect can be calculated 

in the form of linear function A + Bx, which is the local background for a 

given territory. When making quantitative calculations, it is absolutely 
necessary to make allowance for the effects of the local background. After 
allowance has been made for the local background, the calculated anomaly assumed 
another form and exhibited much better agreement with the observed anomaly. 


The last stage of the construction was solution of the inverse problem. 
The elementary projections characterizing contacts between rocks of different 
densities were arranged in a pattern permitting the best possible approximation 
of the anomaly calculated from the selected cross-section to the observed 
anomaly. The results of solution of the inverse problem are shown in Figure 29. 
As we See, within the limits of the possible measurement error, the agreement 
between the anomalies is sufficiently good. However, it can be found immediately 
that in some cases a discrepancy is caused by insufficiency of the geological 
information taken into account in setting up the geological cross-section of 
the first and second approximations. This applies principally to the minor 
details of the cross-section. The large features are characterized quite clearly. 


Even with insignificant initial information about the geological structure 
and with few details about the anomaly as concerns its profile, the system has 
made it possible to obtain good results for all the questions examined. 
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Geological Mapping of Complex Regions 


We will examine a small area within the limits of the Central Bug River 
Valley. A large amount of factual geological and geophysical material was 
collected. The initial data for the selection were the observed gravitational 
force anomaly and the primary geological cross-section containing the density 
characteristics of all the various types of rocks. The calculation was carried 
out by use of a system of parallel profiles intersecting the basic geological 
structures of the region, and also by using small profiles within the limits 
of the individual structures. This distribution of profiles makes it possible 
to correlate spatially both the larger structures and the individual extended 
bodies. In the area under study there are geological features the configurations 
of which have been fairly well studied by geological and geophysical methods. 

In this case the dimensions of the bodies were assigned in advance and could 
not change during subsequent calculations. Individual areas of the Tarnovatskiy 
syncline which has been drilled in detail can serve as an example. The syncline 
is filled with rocks of the metabasite and ultrabasite complex, and these rocks 
differ sharply in their physical properties from the relatively monotonic mass 
of mignatites. The latter are more common in the region in question and serve 
as a kind of intrusive mass for all the other rock complexes. After the direct 
problem had been solved on the basis of the initial geological cross-section, 
the various items of supplemental information and corrections were introduced 
into the model. The amount of work required by this stage depends completely 

on the reliability of the initial geological model. 
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Figure 30 shows a diagram of the geological structure of one of the areas 
in the region being investigated. The Tarnovatskiy syncline is clearly to be 
distinguished. In the northern parts it is a shallow (150-200 meters) fold 
with a very gently sloping of the limbs. To the south one observes subsidence 
of the fold bend to a depth of 600 meters and more, while the width of the 
structure simultaneously decreases. This causes gradual increase in the steep- 
ness of the fold limbs until inverted rock bedding is formed in the eastem 
limb, The structure is made up of dense (o = 3.04 g/cm*) rocks of the metabasite 
complex, but to the south strata of amphiolous gneisses (density 2.9 g/cm*) are 
included in its structure. Small bodies of serpentinites (density 2.46 g/cm) 
are clearly fixed. They cause local minima in the zones of the gravitational 
force maxima. These maxima are caused by the dense rocks of the productive 
stratum, On analyzing the gravitational force field, we found a basis for 
assuming the presence of several small serpentinite bodies in the shallow 
synclinal folds to the west of the Tarnovatskiy syncline. 


It has been established that the bodies of the serpentinites and the 
serpentinous ultrabasites coincide in space with the metabasite massifs and 
normally comprise small concordant massifs of slight thickness. Relatively 
large areas of individual bodies on the map within the limits of the Tarnovatskiy, 
Kapitanovskiy, and Derenyukhinskiy [Translator's Note: as Transliterated] 
synclines are caused by the very gentle slope of the rock strata within the /102 
limits of these structures. 


Thus it has been possible to distinguish within the plan and to observe 
to significant depth a considerable number of large and small folded structures, 
as well as a large number of individual extended bodies which take part in the 
formation of the folds. 


Small geological bodies are distinguished in the region in question. They 
are made up of very dense and highly magnetic rocks of the skarn type which, 
judging from the results of geological studies, are confined to the contact 
zones between rocks of the metabasite and hyperbasite complexes and the granites. 
In rare cases the thickness of these bodies reaches 100 meters over a considerable 
extent. Their formation should probably be associated with the intensive in- 
troduction of free magnetite into the active contact zone. 


Even in the areas where the gravitational force anomaly is little differ- 
entiated above anticlinal structures, it is possible to determine the spatial 
position of individual bodies of dense rock which are included in the structural 
formation. The geological charts and diagrams known up to the present have 
shown large massifs of charnockites in the nuclear portions of these structures, 
the presence of which explains the high background (up to 1000 y) and the mosaic 
character of the magnetic field. There are no charmockites in our massif con- 
struction, and the nature of the magnetic field is explained by the complicated 
structure of these regions, within the limits of which considerable saturation 
is noted by small bodies of dense, intensively magnetic rock such as amphi- 
bolites, pyroxenes, amphibolous and biotitic gneisses, charnockites, and pink 
magnetoactive fine-grained granites. The high magnetic intensity of the rocks 


98 


mentioned is caused primarily by their considerable (from 1 to 3%) content of 
free finely dispersed magnetite. This mineral is observed in microsections. 
The rocks which make up the small bodies within the limits of the anticlines 
are characterized primarily by the presence of free magnetite. These same 
rocks in synclinal structures are manifested as practially non-magnetic ones. 


It is probable that the singling out of free magnetite is caused by the 
higher level of granitization of the rocks of the lower structural level. 
Xenoliths of these rocks have a periclinal closure on relatively sharp (60-70°) 
drop from the center. This has served as the basis for distinguishing these 
structures as anticlines or dome-shaped elevations. In our constructions they 
represent formations in the lower structural stage in the form of hard blocks, 
to which the development of linear folding by the rock in the upper structural 
level is subject. 


In order to map the rocks correctly, one must know the tectonic structure 
of the region under investigation. The fracture tectonics, in causing the 
block-like structure of the region, create a unique regional background for the 
small structures in the uppermost layer of the Earth's crust, and these small 
structures determine the intensive differentiation of the observed gravitational 
force field. The exclusion of or allowance for this background makes it possible /104 
to carry out quantitative interpretation of anomalies properly both for individual 
profiles and for the overall map. In the region under investigation, the border 
of the large blocks, which differ insignificantly in density, are clearly to 
be seen. These tectonic contacts of the fault type are traced to a considerable 
degree in extent and in depth. 


The structural tectonic map of the region is predetermined by the presence 
of two systems of tectonic distortion, a northwestern and a northeastern one, 
of which the first exerts a more significant influence on the formation of 
the small-block nature of the region's precambrian formations. In this process 
a scalariform elevation (of the horst type) is noted in the central part of 
the rigid blocks. It is possible that these elevations cause the formation of 
the gently sloping folds made up of rocks of the gneissic mass. It should be 
noted that not only are the folding of the rocks in the upper structural levels 
subordinate to the contours of the rigid blocks; the tectonic distortions also 
exhibit a tendency to bend around the central portions of these blocks. 


It is quite typical that the elevated blocks consist of lighter rocks. 
These block movements have affected chiefly the structure of the intrusive rock 
mass, which is represented by monotonic migmatites with numerous xenoliths of 
denser rocks. Therefore, it is possible to draw the conclusion that with depth 
the density of the rock of the intrusive rock mass decreases to a certain limit. 
This pattern can be explained by the higher degree of granitization, and accor- 
dingly by the decrease in the number of xenoliths in the dense rock. It is 
probably this that explains the considerable variations in the density of the 
migmatites (2.51 - 2.67 g/cm). 


The examples cited above show the effectiveness of using the minimization 
method when solving problems in structural geology in various situations and 
under various conditions of information coverage. It is natural that the more 
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the area is studied from the geological viewpoint, and the more reliable is the 
first approximation model, the more effective will be the use of an automated 
system. However, even with the most scanty geological information, the system 
makes it possible from the quantatitive point of view to estimate the reliability 
of any assumption and to select the optimum geological model from among many 
optional versions. This latter factor can serve as a starting point for bringing 
geological investigations up to date. 


Conclusion /105 


In this work several questions relation to application of the minimization 
method to interpret gravitational anomalies have been examined. A fairly com- 
plete automated system of interpretation has been obtained: estimation of the 
initial version of the hypothetical geological model, selection of the back- 
ground function, refinement of the density parameters, and refinement of the 
configurations of geological features. Over a period of several years this 
method has been tested with a variety of examples. Good results have always 
been obtained. It may be noted that on occasion, when the geological first 
approximation model was very incomplete, results were obtained which contradicted 
the geological premises. An analysis of this solution has always enabled the 
interpreter to restructure his hypothetical model properly and to overcome any 
contradictions which have occurred. 


The same method can be used to create an automated system for interpreting 
Magnetic anomalies [50]. Of course, the methods described can be used most 
effectively in cases in which the boundaries of the perturbing bodies form 
steep contacts. 


Lastly, the system described can be considered as a subsystem of a more 
general system which would include other calculation methods. 
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Figure 30. Geological Model and Cross-Sections of One 

of the Areas in the Central Bug River Valley: 1) Observed 
Gravitational Force Anomaly, 2) Anomaly Calculated From 
Selected Cross-Section, 3) Anomaly Due to Fracture Tectonics, 
4) Excess Density Value, 5) Serpentinites, 6) Amphibolites, 
7) Amphibolous Gneisses, 8) Biotitic Gneisses, 9) Pyroxenite 
Gneisses, 10) Charnockites, 11) Granites, 12) Migmatites, 

13) Tectonic Distortions. 
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